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Abstract 

The BCS-BEC crossover has received much attention lately, owing especially to its experimental 
realization with trapped ultracold Fermi atoms. Theoretically, the two limiting situations, of paired 
fermions described by BCS theory in weak coupling and of composite bosons undergoing Bose-Einstein 
condensation (BEC) in strong coupling, can be connected with continuity throughout the crossover. This 
evolution encompasses the unitary limit at intermediate values of the coupling, where the scattering length 
for two-fermion scattering diverges. Several quantities have been measured experimentally and calculated 
theoretically in this context over the last several years, with the notable exception of the Josephson and 
related effects. This is in spite of the fact that the Josephson effect is intimately associated with the 
spontaneous breaking of the phase of the complex order parameter which unifies superconductivity and 
superfluidity. 

In the present paper we aim at filling (at least partially) this gap and investigate the evolution 
of the Josephson and related effects throughout the BCS-BEC crossover, by performing a systematic 
numerical solution of the (time-independent) Bogoliubov-de Gennes equations at zero temperature in a 
fully self- consistent fashion. We consider a stationary and uniform current flowing in the presence of a 
three-dimensional barrier with a slab geometry. This extended geometry is specifically required to reach 
the BEC limit of the crossover, where the formation of composite bosons in terms of their fermionic 
constituents requires consideration of wave vectors with components along all three dimensions. In 
addition, we regard the fermionic attraction to extend unmodified over the barrier region, a situation 
that typically applies to ultracold Fermi atoms. The fully self-consistent solution of the Bogoliubov-de 
Gennes equations in such an extended geometry and coupling range represents a non-trivial numerical 
calculation. The numerical strategies and algorithms we have adopted will therefore be described in 
detail, with the aim of easing further independent studies. 

Several results are obtained by the present calculation. The profiles of the magnitude and phase 
of the gap parameter across the barrier are determined under a variety of conditions. We find that 
the Josephson current is considerably enhanced at about unitarity for all barriers we have considered. 
A related enhancement is also found in the contribution to the total current from the Andreev bound 
states, which stem from the depression of the gap profile about the barrier. The Josephson current-phase 
characteristics (relating the total current J to the phase difference 5(j> across the barrier) turn out to 
evolve from the standard J <x sin 8<f> relation to J a cos 5<f)/2, when the height of the barrier is decreased 
at fixed coupling or the coupling is decreased for given barrier. For vanishing barrier height, we find that 
the critical Josephson current approaches the limiting value predicted by the Landau criterion, which is 
determined by either pair-breaking or sound-mode excitations depending on the coupling value. In the 
BCS limit, we reveal the presence of Friedel oscillations in the oscillatory modulations of the gap and 
density profiles. In this limit, we also emphasize the special role played by the Andreev bound state in 
determining the critical Josephson current in the presence of a barrier. Finally, the stability of the two 
branches, out of which the Josephson characteristics are composed, is analyzed by calculating the energy 
required to produce a given spatial profile of the gap parameter. 
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List of symbols: 

as,«F = bosonic (B) and fermionic (F) scattering lengths 
Ep = kp/(2m) — Fermi energy 
£ = total energy of the system 

—g = bare coupling strength of the fermionic attraction 

J = Josephson current 

L = width of a rectangular barrier 

kp = Fermi wave vector 

(k±, k||) = wave- vector components orthogonal and parallel to the barrier 

m 7 mB = fermion and boson (B) mass 

n(x) = profile of the fermionic number density 

no = bulk value of the fermionic number density away from the barrier 
q = wave vector associated with the supercurrent 
s = sound velocity 

V(x) — one-body potential representing the barrier 
Vo = height of the barrier 

{u v ,v v ) — eigen-solution components of the BdG equations 
x = spatial coordinate orthogonal to the barrier 
Z = strength of a Dirac-delta barrier 

S(j) = phase difference accumulated by the gap parameter across the barrier 

A(x) = spatial profile of the (complex) gap parameter 

Ao = bulk magnitude of the gap parameter away from the barrier 

£o = (map)^ 1 = two-fermion binding energy 

ji = fermionic chemical potential in the presence of a supercurrent 
Ho = fermionic chemical potential in the absence of a supercurrent 
Hb = chemical potential of composite bosons 
a = width of a Gaussian barrier 

2<j>{x) — spatial profile of the phase of the gap parameter 
&(x) — condensate wave function solution of the GP equation 
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1 Introduction 



This paper deals with the solution of the (time-independent) Bogoliubov-de Gennes (BdG) equations for a 
system of fermions with a mutual attraction of arbitrary strength. The BdG equations [see Eqs.© below] 
are the analog of the Schrodinger equation for two-component fcrmionic wave functions in the presence of 
a superconducting gap parameter pQ, and generalize to spatially dependent situations the BCS mean-field 
approch for superconducting fermions [2]. Although these equations can, in principle, be utilized at any 
temperature in the superconducting phase, our interest here in extending their applicability to large values 
of the mutual attraction limits, in practice, their use to zero temperature for reasons that will be clarified 
below. 

In particular, we shall be interested in phenomena involving macroscopic quantum coherence which are 
made manifest by the presence of inhomogeneities in the sample. Typically, we will consider a potential 
barrier that partitions the superconducting sample in two halves, such that a Josephson current can be 
established between the two sides of the barrier. The description of the Josephson and related effects will 
then constitute the principal concern of this paper. By our approach, we will be able to explore these effects 
for parameter ranges that could not be explored before, specifically, when the mutual fermionic attraction 
giving rise to superconductivity is made arbitrarily large. From a technical point of view, this program 
requires us to refine to a considerable extent the self-consistent solution of the BdG equations, in such a way 
that one is able to describe the evolution from the weak-coupling limit (which has been considered so far in 
the literature) to the strong-coupling limit, when the original fcrmionic identity of the constituent particles 
is apparently lost in favor of bosonic entities. 

Quite generally, we shall refer to the (dc) Josephson effect as the occurrence of a stationary superfluid 
current J flowing across a barrier of arbitrary shape in an otherwise homogeneous system, such that a well- 
defined relation is established between the value of J and the difference 8<f> accumulated by the phase of the 
superfluid order parameter across the barrier. This current-phase relation is expected to assume different 
forms, depending on the type and shape of the barrier [3J. These different forms can be regarded as defining 
different regimes of the same underlying physical effect. In this respect, the analysis presented in this paper 
will enlarge the range of regimes which can be attained by the Josephson effect, by considering that the 
fermionic attraction can be varied in such a way that the description of the supercurrent flow evolves with 
continuity from a BCS description in terms of fermions and Cooper pairs to a BEC description in terms of 
condensate composite bosons. 

Our emphasis in such an evolution of the original fermionic system when the strength of the mutual 
attraction is increased, relates to what is nowadays known as the BCS-BEC crossover, from a standard 
(BCS) description when the mutual attraction is sufficiently weak, to a description in terms of composite 
bosons which undergo Bose-Eistein condensation (BEC) when the mutual attraction is sufficiently strong. 
The interest in this evolution has been considerably stimulated by the observation that high-temperature 
cuprate superconductors have coherence lengths much smaller than conventional superconductors, and more 
recently by the full experimental realization of the BCS-BEC crossover with ultracold Fermi atoms. 

We will thus begin by giving a short account of this crossover especially in the context of the physics 
of ultracold Fermi atoms, a field which has been very much active over the last several years and offers 
many opportunities to come for exploring previously not accessible physical regimes. We refer the reader to 
Refs.[4, 5] for recent reviews of this field [see also Ref.[B] for a more specific review of a theoretical approach 
to the problem]. In the rest of the paper we shall then proceed, more specifically, to a detailed description 
of the Josephson and related effects throughout the BCS-BEC crossover, thus filling a notable gap in the 
literature. By this calculation, several novel results will be obtained and their physical implications discussed 
at some length. 
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la. The BCS-BEC crossover and the physics of ultracold Fermi atoms 



In conventional superconductors (like Al, Hg, Sn, • • •), electrons with opposite spin pair to form Cooper 
pairs at temperature below the superconducting critical temperature T c . In these superconductors, the 
average size £ pa i r of the pairs is much larger than the average interparticle distance kp 1 (where kp is the 
Fermi wave vector), such that the product &F£pair turns out to be about 10 3 — 10 5 . The Cooper pairs are 
thus largely overlapping and it is not appropriate to regard them as spin-zero bosons (in this case, in fact, it 
is more appropriate to refer to the correlation between fermions with opposite spins at a certain distance). 
This consideration has been central to the BCS theory of superconductivity [2 and has made it possible its 
successful application to conventional superconductors, because the large spatial overlap among Cooper pairs 
justifies a description based on a "mean-field" analysis of the system. In this respect, the superconducting 
transition in conventional superconductors should be regarded as rather exceptional among second-order 
phase transitions, to the extent that a mean-field treatment can accurately account for the experimental 
observations. 

This framework has changed considerably with the discovery of high-temperature cuprate superconduc- 
tors, for which kp £ pa ir is of the order 5—10 only. The BCS assumption of largely overlapping fermion pairs 
is thus not fully realistic for these materials, requiring the BCS theory to be suitably modified. Specifically, 
for these novel superconductors the coupling between fermions seems to be somewhat intermediate between 
situations with highly overlapping Cooper pairs and with composite bosons spatially separated from each 
other. [Here and in the following, the term "composite" bosons relates to the fact that the dissociation tem- 
perature of these entitites is comparable with the critical temperature for condensation, while for "point-like" 
bosons (like 4 He) the two temperature scales differ considerably.] 

These considerations have promoted the development of a theoretical approach which is able to connect 
the BCS limit as described above in terms of Cooper pairs with the Bose-Einstein condensation of composite 
bosons. Since the evolution between these two (BCS and BEC) limits occurs with continuity without the 
occurrence of an intervening phase transition, the phenomenon is referred to as the BCS-BEC crossover. 

There actually exist some theoretical papers [7J |9] which have dealt with this crossover beforehand, 
motivated by the need of producing a rational description of the evolution between the BCS and BEC limits. 
Subsequently with the discovery of high-temperature cuprate superconductors, the interest in the BCS-BEC 
crossover has surged [TUJ HU H21 H31 E] • It is, however, only with the advent of the experiments on ultracold 
trapped atoms that the physics of the BCS-BEC crossover has been promoted to the attention of the scientific 
community at large, being of interest in such apparently different fields like atomic and molecular physics, 
condensed matter physics, nuclear physics, elementary particle physics, and astrophysics. 

On the experimental side, the main novelty introduced by the use of ultracold trapped Fermi atoms 
to explore the BCS-BEC crossover is the possibility of varying essentially at will the attractive interaction 
between fermions of different species [TS], since it is just the occurrence of this attraction that makes it 
possible the formation of Cooper pairs in the medium and of composite bosons in vacuum, starting from 
the isolated fcrmionic species. [We anticipate that the electronic spin quantum number of superconductors 
is now replaced by a different quantum number associated with the atomic hyperfine levels of the Fermi 
atoms.] It is just this possibility that allows ultracold Fermi atoms to be considered as prototype systems in 
Nature with respect to others for which this possibility is strongly limited. 

In practice, in ultracold Fermi atoms the attractive interaction is effectively varied by the use of Fano- 
Feshbach resonances |16j . which are characterized by a resonant coupling between a two-atom scattering 
state with vanishing energy and a bound state in a closed channel. By suitably scanning the position of the 
bound state with respect of this zero-energy level through the variation of a magnetic field, the corresponding 
scattering length of can be varied arbitrarily, passing from negative values before the bound-state is formed 
in the two-body problem, to positive values once the bound state has appeared through the resonance |17j . 

The BCS-BEC crossover has been realized in this way with 6 Li and 40 K ultracold Fermi atoms. The Fano- 
Feshbach resonances utilized for the purpose are sufficiently "broad" that the fermionic many-body problem 
can be effectively described in a simplified manner by a single-channel Hamiltonian with an instantaneous 
interaction [18] . For "narrow" resonances, on the other hand, the single-channel description turns out not 
to be accurate, and it is not even possible in this case to follow adequately the evolution between the BCS 
and BEC regimes with the fermionic densities currently available in the traps [18] . 

Several are the experiments which have been realized so far with ultracold trapped Fermi atoms. Among 
those, we mention: (i) The initial production of composite bosons [TH] [20] ; (ii) The Bose-Einstein conden- 
sation of composite bosons [2T1 [2"2l [2"5] ; (iii) The evidence of a condensate on the BCS side of the crossover 
[2"4l |2"S"1 l2l>] ; (iv) The evidence of a pairing gap [23 [23 [21] ; (v) The measurement of collective modes [301 EI] ! 
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(vi) The evidence for a component in the closed channel [32] ; (vii) The "shot noise" atomic correlations [33] ; 
(viii) The evidence for vortices throughout the whole crossover [34] ; (ix) The introduction of imbalanced spin 
populations [313 [35]; (x) The study of transport properties [37] and band structures [35] in optical lattices; 
(xi) The measurement of thermodynamic properties [35] ; (xii) The measurement of the critical velocity for 
superfluid flow [1^ ; (xiii) The measurement of single-particle spectral features [H] . 

On the theoretical side, an adequate description of the whole crossover bridging the BCS and BEC regimes 
represents a difficult problem which has not been fully solved yet. By its very definition, a single crossover 
theory should be able to capture the two limiting behaviors which aims at connecting (in the present context, 
they are the BCS theory on the one side and the description of a dilute Bose system on the other side.) 
Here, the main difficulty is represented by the need of describing a many-body system in the absence of 
a small parameter, whose presence would allow one to adequately control the theoretical approximations. 
Such a small parameter (identified by the product kp\ap\) actually exists separately in the two BCS and 
BEC regimes, where it can be used to control the approximations to which the given crossover theory 
reduces for dilute Fermi and Bose systems, separately. The lack of a small parameter thus concerns only the 
"intermediate" region between the two regimes (which is centered about the so-called unitary limit where 
| op | diverges) where the product fcp|a_F| can by far exceed unity. A possible strategy in this context consists 
in restricting as much as possible the extension of the intermediate region where a stringent control of the 
approximations is not possible, through a successive refinement of the quality of the approximations to which 
the single crossover theory reduces in the two BCS and BEC regimes that are contiguous to the intermediate 
region. This strategy has been adopted, for instance, by the diagrammatic approach of Refs.[33J 01] which 
includes "pairing fluctuations" beyond mean field, an approach that has produced good numerical results 
also in the intermediate (crossover) region and even at finite temperatures. Note that the inclusion of pairing 
fluctuations is especially required to describe the BEC regime at finite temperatures, for which the density 
of non-condensed bosons acquires a major role. By this approach, fluctuations effects beyond mean field are 
initially considered for a homogeneous system and then extended to the presence of a trapping potential 
that confines the interacting Fermi gas via a local-density approximation, whereby the system is considered 
to be locally homogeneous. 

When dealing with the Josephson and related effects throughout the BCS-BEC crossover of interest 
in the present paper, however, the need arises for dealing with inhomogeneities in the system (such as 
those due to the presence of a barrier separating two superconducting regions) beyond the local-density 
approximation. In this case, the inclusion of fluctuations beyond mean field proves much harder than for 
the corresponding homogeneous case, so that one may be willing to limit to a mean-field approach when 
performing practical calculations. That this may not be too severe a limitation, however, stems from the 
consideration that, in the presence of a spatially dependent external potential, the outcomes of a mean-field 
calculation should already be highly nontrivial, to the extent that the imprint of the excitation spectrum of 
a quantum- mechanical inhomogeneous system can be already found in its ground state [45j . In the following, 
we shall thus limit ourselves to consider the BdG equations at zero temperature and explicitly verify this 
subtle physical expectation from the results of our numerical calculations. 

When calculating physical quantities across the BCS-BEC crossover for the homogeneous case, one may 
often rely on analytic results in the two BCS and BEC limits, against which the results of the fully numerical 
calculations can be confronted. A similar need for confronting with analytic results shows up when solving 
numerically the BdG equations in the presence of inhomogeneities. Special care will thus be devoted in the 
following to identify alternative and independent procedures, which provide independent numerical tests for 
the results of the BdG equations both in the BCS and BEC limits. 

The simplest description of the BCS-BEC crossover occurs at the mean-field level for a homogeneous 
system in the zero-temperature limit. In this case, the usual BCS equation for the gap parameter Ao is 
supplemented by the equation for the density n which specifies the fermionic chemical potential /io, and the 
two coupled equations are solved simultaneously for given density and coupling. These two equations read: 



where k is a wave vector, m the fermion mass, £k = k 2 /(2m) - fi , and E k = + Ag. Note that the 
gap equation ([T]) has been suitably regularized to remedy for the presence of an ultraviolet divergence that 




(1) 




(2) 
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originates from the use of a contact potential to represent the fermionic attractive interaction. We shall 
further comment on this point in sub-section lb below. 

Adopting the Fermi wave vector kp = (S^n) 1 / 3 as the unit of |k|, one sees from the right-hand side of 
Eq.([T]) that the quantity (/cf^f) 1 plays the role of the coupling parameter of the theory. Depending on 
the sign of ap, this parameter ranges from (kp of) 1 ~ — 1 characteristic of the weak-coupling BCS regime 
when aF < 0, to (kp ap)^ 1 ~ +1 characteristic of the strong-coupling BEC regime when ap > 0, across the 
value {kp ap) -1 = at unitarity when \ap\ diverges. In practice, the "crossover region" of most interest is 
limited to the interval —1 ~ (kp ap)^ 1 ~ +1, with the approximate values of the boundaries depending on 
the specific physical quantity at hand. 
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Figure 1: (a) Gap parameter Ao and (b) chemical potential /io for a homogeneous system vs the coupling 
parameter (/cfcif) -1 , evaluated within the mean-field approximation at zero temperature across the BCS- 
BEC crossover. 



The equations (JTJ) and © can be solved analytically for Ao and /io vs (kpap) -1 in terms of elliptic 
integrals, as shown in Ref. [25] • These solutions are plotted in Fig[TJ where we note that a different nor- 
malization has been adopted for /io depending on its sign (namely, the Fermi energy Ep — k F /(2m) when 
fio > and the two-fermion binding energy £q = (map)^ 1 when /io < 0). While Ao is a monotonically 
increasing function of (kp ap) -1 , /io is seen to drop rather sharply from the BCS to the BEC values across 
the crossover region of limited extension. In reverse, one may say that it is just this characteristic behavior 
of the chemical potential to drive the BCS-BEC crossover from the presence of an underlying Fermi surface 
(represented by Ep) to the occurrence of pre-formed fermion pairs (represented by £o)- 

The values of Ao and /io reported in Fig[T]for a homogeneous system will also be of use for the solution 
that we shall undertake of BdG equations in the presence of a barrier of finite width. This is because 
the value Ao will be asymptotically reached far away from the barrier by the magnitude of the spatially 
dependent gap parameter, while the value /io will remain unaffected by the presence of the finite barrier 
which perturbs the density profile only locally. 

lb. Single-channel Hamiltonian and regularization of the fermionic attractive 

interaction 

It is worth commenting more extensively at the outset on the use of a contact potential to represent the 
attractive fermionic interaction, as well on the procedure we have adopted to regularize the divergences that 
may arise in physical quantities when using this kind of potential. 

We have already mentioned that the attractive interaction between two Fermi atoms ( 6 Li and 40 K are 
commonly used in experiments) is provided by a molecular mechanism referred to as Fano-Feshbach resonance 
[16) . This mechanism is characterized by a resonant coupling between a two-atom scattering state with 
vanishing energy in an open channel and a bound (molecular) state in a closed channel. Upon varying the 
position of the bound state with respect to the threshold of the open channel (which is obtained, in practice, 
by varying an applied magnetic field that acts differently on the open and closed channels), the value of the 
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scattering length at the threshold of the open channel can be changed essentially at will, passing from 
negative values in the absence of the bound-state to positive values once the bound state has formed through 
the resonance [T7j. In all cases, an attraction has effectively developed between the two Fermi atoms with 
different internal states (which are conventionally assimilated to "spin up" and "spin down" states). 




-50 1 1 1 ^ — 1 J 

250 500 750 10001250 

B(G) 

Figure 2: Scattering length ap (in atomic units) for the collision of two fcrmionic 6 Li atoms vs the magnetic 
field B (in Gauss) that detunes the scattering from the close channel (the figure is adapted from Ref.[18j). 
The inset show the behavior of aF for the narrow resonance. 



As an example, the scattering length for collisions of two 6 Li atoms is reported in Figj2]vs the magnetic field 
B. Note that in this case ap diverges at the position of the "broad" resonance at 834 G and of the "narrow" 
resonance at 543 G. 

Provided ap is the only relevant parameter which characterizes the low-energy resonant scattering, the 
effective attraction between two fermionic atoms of different spin components at spatial positions r and r' 
can be conveniently represented by a two-body contact potential of the form — g8(r — r') with negative 
strength (g > 0). In this case, some care should be exerted when relating the value of g to the scattering 
length ap which is meant to represent. This relation is obtained by solving the integral equation in the 
low-energy limit for the two-body t-matrix associated with the contact potential [17], which reads: 

to 1 r ko dk to 

^o~f ~ ~~g + J (2^)3 k 2 " ' ( ' 

In this expression, the delta- function potential has been suitably regularized by introducing an ultraviolet 
cutoff fco in the (otherwise divergent) integral on the right-hand side. Equation ([3]) provides a prescription 
to relate the cutoff fco and the strength g to a given value of ap, by letting fco — * oo and g — > in such a way 
that 

9 = — t + r? • 4 

to fco map k{j 

In some circumstances, Eq. (|3J) as it stands can be directly used to regularize expressions with the same kind 
of divergence, as it was already done for the gap equation (TT]). 

With these provisions, one says that the scattering problem for two Fermi atoms, which had originally 
involved (at least) the two channels connected via the Fano-Feshbach resonance, has been effectively modeled 
by a single- channel Hamiltonian with a two-body point-contact potential. The ensuing treatment of the 
many-body Hamiltonian is in a sense "universal" , insofar as it will not depend on microscopic details that 
would distinguish different physical systems [35]. We recall that a similar attitude was adopted when the 
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microscopic BCS theory of superconductivity was formulated in terms of the Gorkov point-contact potential, 
which retains, however, a finite upper cutoff given by the Debye frequency and is therefore limited to dealing 
with weak-coupling situations [2J. 

It was shown in Ref.[18 that ap can indeed be considered as the only relevant parameter characterizing the 
low-energy scattering, whenever the Fano-Feshbach resonance is sufficiently "broad" that other characteristic 
lengths (like the so-called "effective range") can be neglected for all practical purposes. In these cases, it 
turns out that the probability of finding the two-atom system in the closed channel is sufficiently small, such 
that only the open channel can be retained for an effective description of the system. The presence of the 
closed channel remains, of course, essential for an ab initio determination of the scattering length clf itself 
as a function of the magnetic field B, like in the calculation shown in Figf2J In the single-channel model, 
on the other hand, the correspondence between the values of of and B can only be fed in as an external 
information. This is the model that we are going to adopt in the following for a numerical solution of the 
BdG equations across the BCS-BEC crossover in the presence of a barrier. 

lc. Summary of the main results about the Josephson effect throughout the 

BCS-BEC crossover 

When dealing with ultracold Fermi atoms, typical values for the width of the barriers that can be 
introduced in the system are about 5-8 fim and the largest spatial extension of the atomic cloud is about a 
few hundreds /im. Under these conditions, the external magnetic field driving the Fano-Feshbach resonance 
cannot be varied in the interior of the barrier region, in such a way to make the inter-particle attraction to 
vanish there. For these reasons, we shall limit in the following to consider an SsS barrier, with the same 
kind of superfluid extending on the two sides (S) of the barrier and with the fermionic attraction stretching 
unmodified over the barrier (s) region. This differs from the SnS situation encountered while dealing with 
the Josephson effect in more conventional condensed-matter superconductors, whereby a normal (n) region 
lacking the fermionic attraction is embedded in an otherwise infinite superconductor (S). 

An additional property which could be readily considered for ultracold Fermi atoms is the density im- 
balance between the two spin populations [35, 36 . This imbalance acts as an effective magnetic field which 
affects the spin but not the orbital degrees of freedom, and may thus let one to reveal effects which are 
obscured in condensed-matter systems, even as far as the Josephson and related effects are concerned. In 
this paper, we shall limit to consider balanced pupulations only (as it was already done in Eqs.([T]) and (fJJ ), 
and defer consideration of imbalanced populations to a subsequent study. 

This paper will thus be devoted to a systematic self-consistent solution of the time-independent BdG 
equations with an SsS barrier and equal spin populations throughout the BCS-BEC crossover [49]. Focus on 
this crossover will inevitably require us to avoid approximations which are specific of the weak-coupling (BCS) 
regime as usually considered in the BdG literature. In addition, to acquire definite confidence in our numerical 
calculations, we will test them by comparison with calculations based on alternative approaches which are 
specific to the weak-coupling (BCS) and strong-coupling (BEC) regimes. In particular, the connection which 
emerges in the strong-coupling limit between the fermionic BdG equations and the Gross-Pitaevskii (GP) 
equation for composite bosons, that was established at a formal level in Ref.[50j and will be verified here by 
numerical calculations, will enable us to connect phenomena of macroscopic quantum coherence occurring 
in fermionic and bosonic systems, by following their evolution from the one to the other end. Especially in 
this limit, the need for a self-consistent solution of the equations (both BdG and GP) will be evident. 

Besides the spatial profiles of the (magnitude and phase of the) order parameter, our calculations will 
produce the Josephson "current vs phase characteristics" for different couplings and barrier (height and 
width) values. In particular, the value of the current will be analyzed in terms of a complete set of wave 
functions which are solutions of the BdG equations, in such a way that its partition between the contributions 
of the bound and continuum states will naturally emerge from the calculation. Physically, the presence of 
these bound states (called the Andreev-Saint- James states [5T]) is related to the depression of the magnitude 
of the order parameter about the repulsive barrier, and has been strictly associated with the very occurrence 
of the Josephson effect in the literature dealing with weak coupling only [52) . 

Quite generally, we find that the maximum Josephson current which is attainable for given barrier as 
well as the contribution of the Andreev bound states are maximal at (about) unitarity, which thus appears 
the coupling range to be preferred for studying the occurrence of macroscopic quantum phenomena. We 
will attribute this feature to the reduced effect at unitarity of the elementary excitations in the system 
(namely, the pair-breaking excitations on the weak-coupling side and the sound-mode excitations on the 
strong-coupling side), which would act to destroy the superfluid flow. 
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The occurrence of these two types of excitations is explicitly revealed in the limiting situation of vanishing 
barrier height, whereby the infinitesimal barrier acts as a seed for triggering the collapse of the superfluid 
flow according to the Landau criterion [53] . We will show how our self-consistent calculation is able to detect 
the presence of this upper boundary for the maximum value of the Josephson current across the BCS-BEC 
crossover, by signalling the lost of the self-consistent solution when approaching from below the boundaries 
related to the two (pair-breaking and sound mode) types of excitations. These two boundaries, which cross 
each other at (about) unitarity, delimit the region of stability of the stationary Josephson flow in a current 
vs coupling "phase diagram" . The occurrence of these boundaries has been confirmed experimentally by 
a recent study performed at MIT [3D]. In addition to leading to an identification of these boundaries, our 
calculation determines also how the maximum current approaches its Landau limiting value for vanishing 
barrier, for each coupling across the BCS-BEC crossover. 

Several interesting features are further revealed by our calculation away from unitarity. On the weak- 
coupling side, the profiles of the gap parameter and density appear modulated by the occurrence of Friedel 
oscillations of period 2fci?, and the presence of Andreev bound states below the threshold of pair-breaking ex- 
citations relates to the maximum Josephson current in the presence of a finite barrier. On the strong-coupling 
side, on the other hand, we compare the numerical solutions of the BdG equations with the corresponding 
independent solution of the GP equation for the composite bosons, and reveal how characteristic features 
related to the composite nature of the bosons emerge, when the barrier width is sufficiently small with 
respect to the bosonic healing length that the composite bosons cannot be accomodated within the barrier 
region. 

Finally, the stability of the self-consistent solutions obtained by our numerical calculations will be studied 
by calculating the corresponding total energy (both in the presence and in the absence of the barrier) 
associated with the required boundary conditions. This analysis will lead us eventually to relate the onset 
of the Landau criterion (when the supercurrent exceeds a certain threshold) with the energy instability of 
the self-consistent solutions of the BdG equations against the inclusion of fluctuations beyond mean field. 

Additional material is confined to the Appendices, and includes both technical details as well as some 
clarifications about subtle aspects of the Josephson effect. 



2 Implementing the self-consistent solution of the BdG equations 

The Bogoliubov-de Gennes equations embody the BCS mean-field theory of superconductivity for a system of 
fermions with a mutual attraction and subject to an external potential. At zero temperature, solving for these 
equations corresponds to finding a (variational) broken-symmetry ground state of the many-body system in 
an inhomogeneous situation. Although mean-field like, this solution turns, however, out to be highly non- 
trivial, since quite generally in the presence of a non-trivial geometry the imprint of the excitation spectrum 
of a quantum- mechanical system can be already found in its ground state [45] . 

In this Section, we set up the solution of the BdG equations for a specific geometry appropriate to realize 
a Josephson link with a uniform current. As we allow the mutual fermionic attraction to take arbitrary 
values, we shall avoid adopting a number of approximations often used in the literature, which are peculiar 
to the weak-coupling (BCS) limit of this attraction and fail accordingly to account for the formation of 
composite bosons in the opposite strong-coupling (BEC) limit. Since the ensuing numerical solution of 
the BdG equations will unavoidably be quite involved, we shall rely on definite benchmarks with which to 
compare our numerical results in the two limiting (BCS and BEC) situations. These will be identified with 
the delta-like barrier and the Gross-Pitaevskii equation, in the order, to be discussed separately below. To 
avoid overwhelming here the reader with details, some of them will be deferred to the Appendices A and B. 

2a. The Bogoliubov-de Gennes equations 

The fermionic BdG equations read pQ: 

H(r) A(r) \ / u„(r) \ / u„(r) \ 

A(r)* -H(r) J { v v {v) )-^\ v v {v) ) ■ W 
Here, H(r) — — V 2 /(2m) + V(r) — fj, where V(r) is the external potential and /z the chemical potential (we 
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set h = 1 throughout). The local gap parameter A(r) is determined via the self-consistent condition: 



A(r) = 5 £u„(rK(r)* [1 - 2f F (e v )] (6) 

where /p(e) = {e e '^ kBT ' + l) -1 is the Fermi function at temperature T (ks being Boltzmann constant) 
and g the coupling constant discussed in sub-section lb. In addition, the functions {uj,(r), u„(r)} obey the 
orthonormality condition: 

dr [u u (r)*u v ,(r) + v„(r)*v„>{r)] = 5 VV > (7) 



where the Kronecker delta on the right-hand side is readily generalized to the eigenvalues of a continuous 
set. 

The equations are formally analogous to the Schrodinger equation for a two-component wave function, 
in which A(r) plays the role of an off-diagonal potential. We will argue that the fact that this off-diagonal 
potential has to be determined self-consistently is of crucial importance for sustaining the Josephson effect, 
especially on the BEC side of the crossover. Note further that in Eqs.© we have not explicitly included 
a (Hartree-type) diagonal potential (which would also need, in principle, to be determined self-consistently 
1 ) , since it vanishes owing to the regularization we have adopted for the contact potential. 

In the following, we shall restict to consider only the positive eigenvalues e„ of Eq.©, so that in 
Eq.([S]) can be taken to vanish in the zero-temperature limit. One can show, in fact, that positive eigenvalues 
only are required to describe the stationary current flow below the critical current at given coupling. More 
generally, out of the two "branches" of eigenvalues of Eq. ([5]) one can limit to consider only the one that reaches 
large positive energies, with the argument that it suffices to identify a complete set of wave functions. 

Physical quantities of special concern will be the number density and current associated with the solutions 
of Eq.([5]). They are given, respectively, by the expressions: 

n(r) = 2 ]T [Me„) K(r)| 2 + (1 - /fM) k«| 2 ] (8) 



j(r) = — £{/f(£„) [u„(r)*Vu„(r) - (Vu„(r)>„(r)] 

V 

+ (1 - f F {e v )) K(r)V<v(r)* - (V«„(r))«„(r)*]} . (9) 

Again, only positive eigenvalues will be considered for the sums in Eqs.© and ©, so that fpi^u) can be 
set to vanish in these equations. 

The crucial role played by the self-consistent condition © emerges when considering the continuity 
equation associated with the current ([9]). With the help of the BdG equations ([5]) one, in fact, arrives at the 
result: 

V-j(r) = ^ (2/p(e w ) - 1) Q„(r) (10) 

V 

where in the (apparently) "source" term on the right-hand side we have set 

Q„(r) = 4Im{A(rK(r)X(r)} . (11) 

The self-consistent condition ©, however, restores the validity of the continuity equation V ■ j(r) = 0, since 
the right-hand side of Eq. (|10[) vanishes when entering the expression ([B]) therein. 

The strength —g of the local fermionic attraction is suitably eliminated in favor of the fcrmionic scattering 
length Of, by adapting to the inhomogeneous case a procedure introduced for the homogeneous case |43j 
(see also sub-section lb). In that case, one arrived at the following regularized gap equation (that is, for 
which no ultraviolet cutoff is required): 



dk 



(27T) 



3 



, ,. u(k)v(k)* m 



4-7T OF 



(12) 



in the inhomogeneous case one obtains correspondingly for the same coupling: 
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Multiplying both sides of Eq.fT2|) by A(r) and subtracting the resulting expression from Eci. lfTS)) yields 
eventually: 

= ^(l-2/ F y^(rK(r)* 

^ / (1 _ 2Mk)) u{ ^ v{ M* ■ ( 14 ) 

This is the form of the regularized gap equation that we shall implement in the numerical calculations. 
Roughly speaking, its ultraviolet convergence is ensured by the fact that, for a large enough value of their 
energy, the wave functions {u„(r), v„(r)} in the presence of the external potential reduce to those of the 
homogeneous case (with A(r) replacing Ao). The ultraviolet convergence of the expression Q14[l will be 
further discussed in Appendix A. 

We are interested in the solutions of the BdG equations (J3J that are subject to the self-consistent condition 
and sustain a finite value of the current To this end, let's first consider the corresponding solutions 
for the homogeneous case in the absence of an external potential. In this case, one sets 

u tf (r) - « q (k)e i < k+q >- r , v v (t) -» <; q (k)e*( k - q )- r , (15) 

such that A(r) = A q e 2lqr (A q real) and 




(16) 

where A q does not depend on q insofar as e v of Ea. (|16p remains positive (the term q 2 / (2m) therein can, in 
fact, be reabsorbed by expressing /j, = (iq + q 2 /(2m) in terms of the chemical potential /i in the absence of 
the current j = qn/m where n is the uniform density). Correspondingly, the eigenvectors can be expressed 
as follows: 



Uq (k) 2 



(£ + 6c- q ) 2 
(£ + £ k - q ) 2 + A 2 



A 



2 



^ = (fi + & _> + Aj (17) 

where E stands for the right-hand side of Ea. (fT6|) and £k-q = (k — q) 2 /(2m) — «, such that E + £k-q, too, 
does not depend on q. The positive square root is eventually taken of the expressions (|17p . 

The expressions (|15p and (|16[) could have been obtained directly from the requirement of Galilean invari- 
ance, by connecting the solutions of the BdG equations in the lab frame (in which the current has the value 
j = qn/m) and in the frame moving with relative velocity V = q/m (in which the superfluid is at rest). 
Note that, with respect to the moving frame where the current vanishes, in the lab frame the eigenvectors 
([T5|) acquire the phase tp(r) = q r which is linearly varying with r. The presence of a barrier breaks Galilean 
invariance and adds to cp(r) an additional phase <fi(r) (such that <p(r) = q-r + 0(r)) that varies sharply about 
the barrier so as to keep the current uniform. We pass now to discuss in detail the corresponding solutions 
for a particular geometry. 



2b. Solutions for a barrier with slab geometry 

To be specific, we realize a Josephson link with a slab geometry, whereby a potential barrier V{x) is 
embedded in a homogeneous superconductor which extends to infinity on both sides of the barrier. Although 
the profile of the barrier is one-dimensional (along the x direction), the slab is meant to be fully three- 
dimensional since it extends along the two (y and z) orthogonal directions. This is because on the BEC 
side of the crossover the formation of composite bosons in terms of their fermionic constituents requires one 
to include wave vectors with components in all three dimensions. In this respect, we depart from previous 
treatments of the Josephson effect both for fermions [53] and point-like bosons [55J for which the formation 
of composite bosons was not an issue, even though taking into account the two orthogonal directions renders 
the present treatment considerably more involved. For the same reasons, the solution of the BdG equations 
in a fully three-dimensional geometry was also considered in Ref. [56] while studying the evolution of a dark 
soliton throughout the BCS-BEC crossover (albeit for real wave functions only). 
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Figure 3: Step-like representation of V (x), |A(x)|, and 2<f>(x) over a finite mesh of intervals about the position 
of a (Gaussian) barrier centered in x = 0. Here, |A(x)| is normalized to its bulk value Ao and 2<j>(x) has 
been shifted upward for clarity. 



The barrier we shall consider is of the SsS type, whereby the fermionic attraction extends unmodified 
under the barrier region (this is a typical situation encountered with ultracold Fermi atoms). The more 
conventional SnS barrier, with a normal region embedded in an otherwise infinite superconductor, will not 
be considered in the present treatment. Although the potential barrier V(x) can in principle take a generic 
shape, in practice numerical calculations will be performed either with a rectangular barrier V{x) of width 
L and height Vq or with a Gaussian barrier of the form V(x) = Vq e~ x ^ 2<T ' (where Vq > 0). 

Operatively, we proceed to the self-consistent solution of the BdG equations (J5J) by extending to arbitrary 
values of the fermionic attraction the numerical approach to the scattering problem introduced for weak 
coupling in Ref.|57). By this approach, the profiles of V(x), |A(x)|, and <f>(x) are made piecewise constant 
(i.e., step-like) over a number M of intervals, such that 



(Ai real) when X(-i < x < xt with £ = (1, • ■ • , M). [The points {xi\ need not be equally spaced, and one 
can take xq — — oo and xm = +oo.] A typical step-like representation of V(x), |A(a;)|, and 2</>(x) is reported 
in FigH 

The advantage of this representation is that, within each £-th interval, the BdG equations can be solved 
by elementary methods since the problem reduces to that of a homogeneous system with a uniform potential 
level Vi, while the solutions in contiguous intervals are connected by standard continuity conditions. In 
analogy to Eas. (fT5|) . we thus set for the generic solutions of the BdG equations in the i-th. interval: 



for given energy E, where r = [x, y, z), k = {ki, k y , k z ), and q = (q, 0, 0). In this interval, the BdG equations 
then reduce to: 



V(x) = Vi and A(a?) = A e e 2 ^ qx+M 



(18) 



u(t;q,E) 
v(r;q,E) 



^( g ,^) e j(k+q) - r+ ^ 



(19) 




Eu e (q 7 E) 




(20) 



where we have introduced the notation 



fit = H — Vi 




(21) 



2m 



12 



which identifies the (local) reduced chemical potential. The solution for E in Eos. ([20]) is thus formally 
analogous to the expression (|16p with the local replacements k x — > kg, fi — > fig , and A q — > A^: 




E=^ + \ £r + £r-fc + *2- (22) 



The reason for replacing the variable k x in Eq. (|16p by the local value kg is now apparent. For given values 
of fig and in the £-th interval, there are in fact four distinct (complex) values of kg that satisfy Eq. (j2"2")) 
for fixed E. In the absence of current (i.e., for q = 0), they correspond to the electron-like and hole-like 
excitations discussed in Ref.[58 . The evolution of these values of kg in the complex /c-plane when E varies 
from to +oo (depending also on the sign of fig) will be discussed in detail below. 

Correspondingly, the eigenvectors of Eqs.(|2"0)) can be cast in a form similar to the solutions (fT7|) . Intro- 
ducing the notation £,k t -q = {kg — q) 2 /(2m) — jig, we arrive at the expressions: 

Ug(q,E) 2 ( E + ^) 2 



(E+^ q r + Aj 

V ^ E)2 = 7FH A La 2 (23) 

(E + ikt-q) + 

where again the positive square root has to be taken. Note that we have conventionally imposed on the 
solutions (|2"3")1 of the homogeneous system (|2"U)) the condition ug(q,E) 2 + vg(q,E) 2 = 1, even though this is 
not directly related to the overall normalization condition Q for the total wave function which extends over 
the whole set of intervals. 

Quite generally, to each value of [with n = (1, ■ • • ,4)] that satisfies Eq.(J22), we can associate the 
two-component wave function [x; q, E) 

for xg-i < x < xg. We use the following convention: (i) n — 1 represents an electron-like excitation 
impinging from the left (of the barrier); (ii) n = 2 an electron-like excitation impinging from the right; (iii) 
n = 3a hole- like excitation impinging from the right; (iv) and n = 4 a hole- like excitation impinging from 
the left (the inverted convention for the hole excitations corresponds to the inverted current flow for hole 
with respect to electron excitations). [There will actually occur an exception to the correspondence of this 
convention with the actual physical situation, in the energy region when < for £ = 1,M.] Note that 
the terminology "electron" and "hole" is mantained for historical reasons, even though the fermions here 
considered are ultracold atoms. 

In terms of the wave functions (|2"4"1) , we implement the solutions of the BdG equations in the ^-th interval 
as follows: 

* t (x\ q, E) = ojffi; q, E) + b e T {2) (x; q, E) + cftf } (x; q, E) + d e T ( f\x; q, E) (25) 

which generalize the solutions (|19[) when all possible values of kg are considered, apart from the common 
factor exp(k y y + k z z) which has been dropped out. The coefficients (ag, bg, eg, dg) are determined by applying 
the continuity conditions to the wave functions (|25|) and their derivatives at the points {xg ; £ = 1, ■ • • , ill— 1} 
between the £-th and (£ + l)-th intervals, plus imposing appropriate boundary conditions in the 1-st and 
M-th intervals at the edges. 

The continuity conditions at the point xg separating the £-th and (£ + l)-th intervals can be expressed in 
the compact form: 

Mg(x = xg) Wg = M l+1 (x = xg) W l+1 (26) 



where Wg is the column vector formed by the coefficients of Eq. (|25)l 

ir, = " (27) 



eg 
\dg J 
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and M^(x) is the 4x4 matrix composed by the following column vectors (n = 1, 2, 3, 4): 



/ M t (x)i n \ 
Mg(x) 2n 
M((x) 3n 

\ Mg(x) in J 



i(k^ n) + q)uf\q,E)e i ^ x+ ^ 
V i(k^ n) - q)v^ n) (q,E)e- i ^+M J 



(28) 



(the dependence of Mg(x) and Wg on (q,E) is understood throughout). Since there are (M — 1) continuity 
conditions of the type (|26p (each of them corresponding to four conditions owing to Eas.(|27|) and ((28]) ') and 
AM coefficients ({ag, bg, eg, dg}; £ — 1, • • • , M), the problem is apparently under-determined. As we will see 
shorthly below, however, it turns out that in the 1-st interval only the coefficients (&i, c\) and in the last M- 
th interval only the coefficients (omt^m) need be determined, thus reducing the total number of coefficients 
to 4(M — 1). Grouping then together all continuity conditions in a single equation, one is left with solving 
the equation 

AW = B (29) 
where W is now the column vector with 4(M — 1) components 

/ bi \ 

Cl 



W 

(If; 

a M 
V d M J 

and A is the 4(M — 1) x 4(M — 1) matrix represented in the following block form: 



(30) 



Ml(il) 



V 



-M 2 (ii) 



-M 3 (x 2 ) 
M3( x 3) 



-M 4 (« 3 ) 



M M-2( I M-2) 



- M M-l( x M-2) 
M M _ 1 (as M _ 1 ) 



-M M (x M 



In the above expression, Mi and Mm arc 4x2 matrices constructed from the pairs of column vectors (|28|) 
with n = (1, 4) and n — (2, 3), respectively. 

The vector B on the right-hand side of Eq. (f2"9"]l is specified by the boundary conditions in the 1-st and 
M-th intervals, and takes accordingly different forms depending on the energy ranges identified by the 
expression (|2"2"|) whereby one sets A^ equal to the bulk value Ao, Ve = 0, and kg = k. The ensuing expression 
assumes different shapes depending on the sign of p, = fi — {ky + fcf)/(2m). In Fig|4]we plot this expression 
vs k (with a value of q > for which it remains positive) for given values of Ao and fx. The two cases with 
ft, > [FigfJJa)] and /2 < [FigHJb)] are considered separately. Note that negative values of ft result for 
large enough values of (k y ,k z ) even when [i > 0. 



Four energy ranges are identified when /2 > and three additional ones when fl < 0, in the following 
way: 
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Figure 4: Dispersion relation (|2"2"]l for I — 1 or £ — M (whereby = A<j, Vi — 0, and kt — k), with (a) 
p, > and (b) fi < 0. These figures identify the seven energy ranges for which the boundary conditions in 
the outermost (left and right) spatial intervals are specified. 
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Range 7 
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no 
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Table 1: Allowed form of the wave function (|25p in the outermost left and right intervals (where the boundary 
conditions are enforced) depending on the value of the energy E. 



With reference to the two different shapes of the dispersion relation reported in FigHJ E^ n stands for the 
absolute minimum on the left side of Figflja), -E^n* f° r the local minimum on the right side of FigfJJa), 
B max for the local maximum of FigHfa), E m - ln for the absolute minimum of FigHJb), and Eq for the value 
corresponding to k — in Fig|4jb). 



In Figs E] and [5] we show how the different solutions of Eg. ([!?!?]) move in the complex fc-plane by 
increasing E in each of the above seven energy ranges. From the form (|25p of the wave function in the 
outermost left (£ — 1) and right (£ = M) intervals, we conclude that the allowed (i.e., non divergent) 
components (|2"4")l must correspond to Im{fc(™)} < when I = 1 and lm{k^} > when I — M. This leads 
us to exclude a priori some of the components (|24[) by setting to zero the corresponding coefficients in the 
wave function (|25[) . as conveniently summarized in Table 1. 

Notice that, by allowing k^ to span the complex fc-plane when the energy E varies from to +oo, we 
avoid relying on the so-called Andreev approximation which is typical of a BCS treatment, whereby only 
wave vectors close to the minima of FigQJa) are effectively retained (see, e.g., Ref.|59j). 

To proceed further, we recall that our purpose here is to identify a complete set of energy eigenstates 
that enter the expressions ©, ©i an d © of the physical quantities of interest. Apart from the presence of 
bound states which are localized about the barrier and exist only in the energy ranges 1 and 5, in the other 
energy ranges we are going to use "outgoing boundary conditions''' for waves impinging on the barrier from 
the left and right, in order to identify the contribution to the complete set from energies in the continuum. 
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Figure 5: Paths followed in the complex £;-plane by the four solutions kS n > by increasing E for jl > 0. The 
panels correspond to ranges 1-4 of Table 1. 



We are now in a position to specify the vector B of Eq. ([2T))) in each of the above seven energy ranges, 
which we treat separately as follows: 

(i) According to the first panel of Fig0 in range 1 there is no real solution to Eq.([22]). From Table 1 the 
wave functions in the outermost intervals are then given by: 



*i(a;; q, E) = bxTf\x; q, E) + Cl T^{x; q, E) 



(3), 



(31) 



for Xq < x < X\, and 

*m(x; q, E) = a M r$ (x; q, E) + d M T$ (x; q, E) (32) 

for xm-i < x < xm- In this case, the vector B is identically zero and Ea. ((29l) reduces to a homogeneous 
systems for the coefficients W , which admits nontrivial solutions only when the determinant of the matrix 
A vanishes. This occurs only for particular values of the energy E in this range, which thus identify the 
bound-state energies. [Note that there exists a continuous branch of bound states associated with different 
values of {k y , k z ).] 

(ii) In range 2, a\ — cm — according to Table 1, while according to the second panel of Fig[5] only 
the components (I24[) corresponding to and fc^ 4 -* propagate undamped. The use of outgoing boundary 
conditions then results in two independent types of external wave functions, namely: 
- Hole-like excitation impinging from the left, such that in the outermost intervals 
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Figure 6: Paths followed in the complex fc-plane by the four solutions k^ n > by increasing E for jl < 0. The 
panels correspond to ranges 5-7 of Table 1. 



ty M (x;q,E) = a M T^(x;q,E) + d M T^(x;q,E) 

-(4), 



(33) 



where the coefficient d\ of (x; q, E 1 ) has conventionally been set equal to unity (the overall normalization 
of the total wave function will be determined at a second stage). According to the notation and 
the convention (J2SJ), in this case B — (— Mi (£1)14, — Mi (2:1)24, — Mi (2:1)34, — Mi (2:1)44, 0, • • - ,0) T (where T 
stands for transpose). 

- Electron-like excitation impinging from the right, such that 



®i(x;q,E) 
^ M (x;q,E) 



b l ?? ) (x;q,E) + c l ??>(x;q,E) 
aM^M( x '^^ E ) + t m( 2: ;9: e ) +4fT^(i;g, E) 



(3), 



(34) 



( 2) 

where now the coefficient 6m of T\^{x]q,E) equals unity. In this case, B — (0, • • • , 0, Mm(#m-i)i2, 
M m {xm-i)22, M M {x M -i)32, M M {x M -i)i2) T ■ 

(iii) In range 3, all components (j24]l propagate undamped [cf. the third panel of Figj5] and the four 
coefficients (01, di, 6m j c m) are in principle all different from zero [cf. Table 1]. These coefficients are fixed 
by the outgoing boundary conditions that we have chosen to determine a complete set, which select four 
independent types of external wave functions as follows: 

- Electron-like excitation impinging from the left, such that (ai = 1, d% = 0, 6m = 0, cm = 0), 



Vi(x;q,E) 
V M (x;q,E) 



T\ l \x;q,E)+b 1 rY>(x-q,E)+c 1 r\*>(x-q,E) 



-(2)/ 



-(3)/ 



(35) 
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and B = (-Afi(a;i)ii, -Afi(xi) 2 i, -Mi(xi) S i, -Afi(xi) 41 ,0, • • • , 0) T . 

- Electron-like excitation impinging from the right, such that {a\ = 0, d\ = 0, &m = 1, cm = 0), 

*i(a:;«/,£0 = hT^\x; q, E) + ciTtf \x; q, E) (36) 
^M^;?^) = a M r[ 1 I \x;q,E) + r M ) (x;q,E)+d M r$(x;q,E) , 

and -B = (0, • • • , 0, M m (xm-i)i2, Mm{xm-\)v2, M m (^-1)32, Mm{xm-i)a2} t ■ 

- Hole-like excitation impinging from the left, such that (a\ = 0, d\ = 1, 6m = 0, cm = 0), 

*i(x;q,E) = b 1 r^(x;q,E) + c 1 r^\x;q,E)+T^(x;q,E) 

* M (x;q,E) = a M r^(x;q,E)+d M r^(x;q,E), (37) 
and B = (-Ml (2:1)14, -Mi(^i) 2 4, -Mi (3:1)34, -M x (2:1)44, 0, • • • , 0) T . 

- Hole-like excitation impinging from the right, such that (01 = 0, d\ = 0, 6m = 0, cm = 1), 

*i(x;q,E) = b 1 r ( ?\x;q,E)+c 1 r ( ?\x;q,E) (38) 

and B = (0, • • • , 0, M M (x M --i)i3, Mm (2^-1)23, M M (x M - 1)33, M M {x M - 1)43 ) T - 

In each of the above four cases, one may identify the different processes as being alternatively the normal 
reflection, normal transmission, Andreev reflection, and inverted transmission. [The last two processes 
are characteristic of the superconducting state [51j . and require one to consider the hole- like wave vectors 
(k^ 3 \ k^) besides the usual electron-like wave vectors (fc^\ k^) characteristic of the normal state.] 

(iv) In range 4, only the components ((24)) corresponding to feW and k^ propagate undamped [cf. the forth 
panel of FigJS] and c?i = cm = [cf. Table 1]. It follows that the only allowed excitations are electron- like 
impinging either from the left [with (eti = 1, 6m = 0), cf. Eq. (f3"5)) ] or from the right [with (ai = 0, 6m = 1), 
cf. Eg. (f3"6")) ] . whereby the vector B has the same form as in those cases. 

(v) The situation in range 5 is equivalent to that in range 1 [cf. the first panel of Fig[6]and Table 1]. There 
is no real solution to Eg. (|2"2")1 , the wave functions in the outermost intervals are again given by Egs. (|3"Tj) 
and (|32p. the vector B vanishes identically, and one is looking for bound-state solutions localized about the 
barrier for particular values of E. 

(vi) The situation in range 6 is similar to that in range 2 [cf. the second panel of Figj6]and Table 1]. This is 
because, when the wave vector fc^ 1 ^ is negative, it should be associated with a hole-like excitation impinging 
from the left rather than with an electron-like excitation impinging from the right. The corresponding 
wave function in the outermost intervals is given by the expression (1331) . with in the upper equation 

now replaced by , and B — (— Mi(xi)n, —M\ (2:1)21, —Mi (#1)31, —Mi (3:1)41, 0, • • • , 0) T . The excitation 
associated with k^ 2 \ on the other hand, remains an electron-like excitation impinging from the right as in 
Eq.®. 

(vii) The situation in range 7 is fully equivalent to that in range 4 [cf. the third panel of FigE] and Table 1] . 
Here, the excitation associated with fcW > represents an electron-like excitation impinging from the left. 

There remains to specify how the ortho-normalization condition |[7J) for the above wave functions is 
implemented in the different energy ranges. In ranges 1 and 5 where the wave functions are normalizable 
(being exponentially localized about the barrier), one breaks up the x-integration into M intervals and uses 
the generic form (|25p of the wave function in the internal intervals xe < x < X( + i with I = (1, ■ • • , M — 2) 
as well as the forms (|3"Tj) and (f3"2")> in the outermost intervals. In all other ranges where the energy E is 
continuous, on the other hand, the normalization per unit range of the wave vector k x is obtained by replacing 
the unit coefficient of the impinging excitation by l/\/27r and rescaling the other coefficients accordingly. 
We have also explicitly verified that degenerate wave functions with given E of the kinds discussed above are 
orthogonal to each other, as required when constructing a complete set. Finally, all types of wave functions 
are eventually normalized per unit range of the wave vectors k y and k z by multiplying them by l/(27r). 

2c. Weak-coupling (BCS) limit and the delta-like barrier 

The general computational scheme, which we have set up in detail in sub-section 2b to solve the BdG 
equations, holds quite generally for any coupling across the BCS-BEC crossover from the weak-coupling 
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(BCS) to the strong-coupling (BEC) limits. Implementation of this computational scheme rests unavoidably 
on sophisticated numerical procedures (as described in Section 3) and involves heavy numerical calculations, 
which require that the outcomes of the calculations need be judiciously assessed before extracting firm 
physical conclusions from them. To this end, it is of definite importance to rely on independent benchmarks, 
with which the outcomes of the calculations could be confronted at least in some relevant limits. The 
crossover problem itself identifies these benchmarks on physical grounds, separately for the two limiting 
(BCS and BEC) situations. We discuss here the BCS limit and defer the discussion of the BEC limit to 
sub-section 2d. 




Figure 7: Pair coherence length £ pa ; r (dashed line) and phase coherence (healing) lenght £ p hasc (full line) 
evaluated at zero temperature within the mean-field approximation according to their definitions as given 
respectively in Refs.[T2] and [50], vs the coupling parameter (^af)" 1 across the BCS-BEC crossover. Both 
lengths are in units of the average interparticle distance kp 1 . Note that in weak coupling the two lenghts 
differ by an irrelevant numerical factor (£ pa i r — (3/v2)£phase) owing to their independent definitions, so that 
in the plot the two curves result parallel to each other in this coupling regime. 



We begin by considering the homogeneous situation, whereby two length scales can be identified as being 
relevant to the BCS-BEC crossover problem at zero temperature [SD]. The first one is the pair coherence 
length £ P air for two-fermion correlation, which corresponds to the size of the Cooper pairs in the BCS limit 
and to the radius of the bound-electron pairs in the BEC limit. The other one is the phase coherence lenght 
£phase) which is associated with the spatial fluctuations of the order parameter (and is often referred to as the 
healing length in the bosonic literature). These two lengths coincide with each other (apart from a trivial 
numerical factor originating from their different definitions) in the BCS limit, where they also reduce to the 
Pippard coherence length characteristic of the superconductor electrodynamics [61) . In the BEC limit, on 
the other hand, the two lengths differ considerably, representing the "intra-" and "inter-" boson correlation, 
respectively. 

The two lengths are plotted in Figj7] versus the coupling parameter [UfCLf)^ 1 spanning the BCS-BEC 
crossover. One sees from this figure that in the BCS limit (/c^of) -1 *C —1 both lengths grow without 
bound, exceeding eventually any other finite length scale that may enter the problem. They will exceed, in 
particular, the characteristic width of a barrier (like that considered in sub-section 2b), so that the specific 
value of this width becomes irrelevant when calculating the scattering properties. 

Under these circumstances, one expects on physical grounds that the barrier could be assimilated to a 
Dirac-delta barrier with potential V(x) = ZS(x), where the constant 

/+oo 
dxV{x) (39) 
-oc 

represents the total area of the potential. For instance, Z — VqL for a rectangular barrier of width L and 
height Vq, and Z — Vo<7V2tt for a Gaussian barrier with V(x) — Vq e~ x ^ 2a '. In practice, this expectation 
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implies that alternative numerical calculations with different potential V(x) that share the same area ([3T))) 
should produce the same results. In the following, we shall use this kind of universality as a stringent test 
on the correctness of the numerical calculations in the BCS limit of the crossover. 
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Figure 8: Self-consistent profiles of |A(x)| and 2<fi(x) in the weak-coupling limit ((kpap)^ 1 = —1.5) for two 
rectangular barriers with the same value of Z. Full lines corresponds to (Lkp — 0.01, Vo/Ep = 30.0) and 
dashed lines to (Lkp = 0.5, Vo/Ep — 0.6). The function 2</>(x) has been shifted upward by 0.4 for clarity, 
and the supercurrent flows with q/kp = 0.045. 



As an example, we report in FigJ5]the profiles of the magnitude |A(x)| (in units of its bulk value Ao) and 
phase 2<j>{x) obtained by the self-consistent solution of the BdG equations for the coupling value (fcj?af) _1 = 
— 1.5 and for two rectangular barriers of width L and height Vq, in such a way that the product Z — LVq is 
kept fixed. Although the full line corresponds to a barrier more closely resembling a delta function (having 
L/£phase — 2.9 x 10~ 3 ), one sees from the figure that there are actually only minor differences with respect 
to the dashed line corresponding to the larger barrier (having L/£ p hasc — 0.15), thus implying that the limit 
of a delta function is readily approached as soon as £/£ p hasc <C 1. In Section 4 we shall return to this issue 
of how our numerical calculations approach this universal behavior. 

An additional interesting issue, which has been discussed over the time in the literature when dealing 
with the BCS (weak-coupling) case, is the role played by self-consistency when solving the BdG equations in 
this limit. Contrary to the BEC (strong-coupling) limit for which self-consistency always needs to be fully 
implemented, in the weak-coupling limit there are, in fact, situations for which the self-consistent process may 
not be necessarily activated in the barrier region. These cases, when the solution of the scattering problem 
can be found analytically without requiring self-consistency, thus represent additional valuable tests for the 
fully numerical calculations. 

Restricting to a Dirac-delta barrier as above, the non-self-consistent calculation is conveniently defined 
by taking |A(x)| = Ao everywhere (except, possibly, at x = where we may set A(x) = 0), while 4>{x) = 
for x < and 2</>(x) = 8<j) for x > 0. In the non-self-consistent calculation this ansatz is used to solve the 
BdG equations only once, without recalculating the profiles of |A(x)| and 4>(x) from the solutions of these 
equations. In such a case, it can be shown that only the bound-state solutions contribute to the current 62J. 
These solutions have a form similar to Eqs. (|3"Tj) and (j3"2"|) , namely, 

= 6l ( J ) e"« W - + c L ( J ) e* W - (40) 

for x < 0, and 

= *R y - e - lSH 2 )z %k + d R y _ e _^ /2 je* ■ (41) 
for x > 0, where k^ = J2m{fx + V A q - E 2 ) and k^ = W / 2m(/I - iy/ 'Ag - E 2 ) (with jj, = /i-kjj /(2m) > 
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and E < A ), while 



\ 



\ 



1 i-i^g 

2 I £ 



(42) 



(43) 



By comparing with Eqs. ([3"T]) and ([3"2"|) . note that in the non-self-consistent solutions that we are consid- 
ering we have set q = 0, so that the (complex) wave vectors (k^\ k^ 2 \ k^ 3 \ k^ ) reduce in the order to 
(fe( e ),-fe( e ),A;W,-fe( h )) as defined above. 

In the weak-coupling limit, Ao <C fi — Ep where Ep — k F /(2m) is the Fermi energy. We also anticipate 
that E ki Ao for the bound-state solutions we are looking for, so that the real parts of k^ and fcW are 
approximately equal to k^ — y/2mjl ~ \J kp — k| in the weak-coupling limit we are considering, while their 

imaginary parts are much smaller than the real parts and will thus be neglected in the following calculation. 
This corresponds to the so-called Andreev approximation, which has invariably been adopted in the solution 
of the BdG equation in weak coupling (see, e.g. Ref.[6"3"]). 

Applying the appropriate boundary conditions at x = 0, namely, \& p{x = 0~) = ^ r(x = + ) = ^(x = 0) 
and d 1 ^ / dx\ - — d^f r(x) / dx\ + = 2mZ^(x = 0), we arrive at the following homogeneous system for the 
unknown coefficients (&l, cl, an, da): 
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(44) 



where we have set Z^ — mZ/kp. It can be verified that the only solution to Eq. (|4"4"]) has the form: 



E ( 



Ao 



COS( 



\ 2(1 + % 



(45) 



This expression can be further manipulated by introducing the transmission coefficient Tk~ for the normal 
case, which equals (1 + 2?) _1 for a Dirac-delta barrier. One obtains eventually: 



E ( 



kp) = Aoa/1 - T fe , sin 2 ,50/2 



(46) 



Under these circumstances, the Josephson (current vs 5(f>) characteristic can be obtained [53] by taking the 
derivative of the above expression: 



J(S4>;kp) = -2 



dE Q ( 



A T fe 



1 - T kp sin 2 S0/2 



(47) 



Note that the energy eigenvalue (|46[) and the associated current f4"Tj) depend on the wave vector k|| parallel 
to the surface of the barrier through kp,, so that the expression ([4"T]) has to be integrated over k|| to produce 
the total current. This integration will be done in a closed form in Appendix B. 

The calculation simplifies in the limit of strong (mZ /kp 3> 1) and weak (mZ/kp <C 1) barriers. For a 
strong barrier, TV — Z^ 2 <C 1 and 



</(&/>; kp) 
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2i? 



2m 2 Z 2 



(48) 



such that 



J{54>) = 
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2^ 



c?fc|| fc|| J(S(f>; kp,) 
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(49) 
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where we have introduced the dimensionless quantity Zp = mZ/kp. We recognize in this type of J vs 
dependence the standard form of the dc- Josephson current [65] . 

This form can also be reproduced numerically from the full self-consistent calculation, by considering 
a sequence of rectangular barriers of width L (such that L/£ p hase <C 1) and height Vq with progressively 
increasing values of Zp = (Vq/ Ep)(Lkp)/2 such that Zp ^> 1 eventually, and fixing to the value 50 of 
Eq. p5)l the asymptotic phase difference 2 [0(x — +oo) — 0(x — — oo)] accumulated by the gap parameter 
across the barrier. 
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Figure 9: Self-consistent profiles of (a) |A(x)| and (b) 20(x), and (c) Josephson characteristics J(S4>) for the 
coupling (kpap)' 1 — —1.5, obtained for a rectangular barrier with Zp = 0.25 (full lines) and 0.5 (dashed 
lines). In the three panels, the dotted line represents for comparison the result of the non-self-consistent 
calculation. The asterisks in the lower curve of panel (c) further represent the results of the self-consistent 
calculation with Zp = 2.5 for which the limit 1 <C Zp has effectively been reached, so that no deviation can 
be appreciated from the non-self-consistent calculation (dotted line). 



A comparison between self-consistent and non-self-consistent calculations is reported in Fig[5J We see 
from this figure that in the limit Zp 1, not only the Josephson characteristics J(0) obtained by the 
fully self-consistent calculation approach the result (|49p obtained analytically by the non-self-consistent 
approximation, but also the profiles of |A(x)| and 20(x) approach the corresponding non-self-consistent 
profiles (barring, of course, the details manifesting at the microscopic length scale kp 1 that are associated 
with the quantum nature of the fermion gas). The reason is that, when Zp ^> 1, the self-consistent profile of 
|A(x)| is strongly suppressed locally in the barrier region over a length scale of the order of Cy — (2mVo) _1/ ' 2 
or kp 1 (depending which is the largest one), whereas outside the barrier |A(x)| approaches its asymptotic 
value A over the larger length scale £ p hasc- At the same time, when |A(x)| gets close to zero for x = 0, 
the phase 4>(x) is allowed to vary sharply about x = without making the complex number \A(x)\e l2 ^ x ^ to 
move appreciably in the complex plane. The result (|49[) is thus universal in the weak-coupling limit, since it 
is obtained when both conditions L/£ p hasc "C 1 and £y /£phase "C 1 are satisfied, leaving £ p hasc as the only 
relevant length scale in the problem (again, apart from the microscopic length scale kp 1 ). 

Returning to the general expression (|47p. we now consider its limit for a weak barrier when Zp <C 1. In 
this case we approximate Tk~ ~ 1 in Eq.(|4"7[). such that J(50; kp) ~ Ao sin 50/2 and 

J{W) = 7~~ / F dk« fcii J(8<t>;kn) ~ %^ sin<W2 . (50) 
2tt J An 

This result can be obtained more thoroughly (as shown in Appendix B), by first integrating over ku the 
generic expression (jT7|) which is valid for any barrier strength, and then taking the limit of the result for a 
weak barrier. Contrary to the case of a strong barrier discussed before, however, for which we have shown 
that the self-consistent and non-self-consistent calculations lead essentially to the same results, in the case of a 
weak barrier the non-self-consistent result (j50|) is completely unreliable because it has the wrong dependence 
on 6<j) and too large a prefactor as well. It will, in fact, be shown in Section 4 that for a generic weak barrier 
(which is not necessarily of the Dirac-delta type) the Josephson characteristics are proportional to cos 50/2 
and not to sin 50/2, and further that their prefactor cannot exceed the one dictated by the Landau criterion, 
which equals Aokp / (3ir 2 ) in the weak-coupling limit we are considering in this Section. 
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2d. Strong-coupling (BEC) limit and the Gross-Pitaevskii equation 



In the previous sub-section we have emphasized the importance of relying on definite benchmarks, with 
which the outcomes of the numerical calculations can be confronted in some relevant limits at least. We have 
also pointed out that for the BCS-BEC crossover these benchmarks have to be looked for in the two limiting 
BCS and BEC situations. In particular, in the BCS (weak-coupling) limit we have argued that the results of 
the fully numerical calculation are bound to depend only on the total area of the potential associated with 
the barrier and not on its specific shape, and we have consequently identified the reaching of this kind of 
universality (which is internal to the numerical solution of the BdG equations) with the desired benchmark 
in this limit. 

In the BEC (strong-coupling) limit that we pass now to discuss, it turns out that one can do even better 
than this and identify a benchmark which is independent from the numerical solution of the BdG equations 
in the presence of a potential barrier. We may again start our reasoning from the homogeneous situation 
where, as discussed in the Introduction, a description in terms of composite bosons naturally emerges in 
the strong-coupling limit of a fully fermionic (zero-temperature) mean-field approach provided the chemical 
potential is suitably renormalized. In the presence of inhomogeneities (such as a spatially varying external 
potential) which induce spatial variations also in the gap parameter, the fermionic approach is embodied by 
the BdG equations |J5J) that generalize the homogeneous BCS approach. In the strong-coupling limit of these 
equations, we still expect on physical grounds a description in terms of composite bosons to emerge even in 
the presence of an external potential, possibly with some limitations on its spatial variation. In addition, to 
the extent the internal structure of the composite bosons is irrelevant that they can be treated as point- like 
bosons, one expects these bosons to be described by the Gross-Pitaevskii equation [15] [55J in the presence 
of a spatially varying external potential acting on the bosons. 

In this context, it was shown in Ref . [50] that, in the strong-coupling limit 1 <C (fcp-ap-) -1 of the fermionic 
attraction, the fermionic BdG equations fl5J can be mapped onto the GP equation 

_Zl$( r )+2T/(r)$(r) + ^^|<i>(r)| 2 a>(r) = » B $(r) (51) 

for the condensate wave function <E>(r) of composite bosons with mass 2m. In the above equation, V(r) 
is the same potential of Eqs.([5]) and \ib = 2/i + Eq is the chemical potential for composite bosons (where 
£o = (map) -1 is the two-fermion binding energy, with fis -C £o)- The two functions $(r) of Ea. (j5Tj) and 
A(r) of Eq.Q are related via the expression 

*(r) = A(r) sj^f- ■ (52) 

We emphasize that this mapping is established at low-enough temperatures, so that all composite bosons are 
condensed and the mean-field approach is appropriate to describe (at least qualitatively) the whole BCS-BEC 
crossover. 

The derivation of the GP equation ((STj) from the BdG equations ([5]) was achieved in Ref. 50 , by trans- 
forming first the differential BdG equations into the associated coupled integral equations for the normal 
and anomalous single-particle Green's functions, and by expanding then the latter equations in terms of the 
small parameter A/|/i|, in such a way that only integrals containing the non-interacting Green's function 
subject to the same external potential remain to be evaluated. In this process, a crucial step was represented 
by the introduction of a local-density approximation for the non-interacting Green's function, which plays an 
analogous role to the eikonal approximation introduced by Gorkov to derive the Ginzburg-Landau equation 
near the critical temperature in the opposite BCS (weak-coupling) limit [57], whereby the presence of an 
external magnetic field is responsible for the spatial inhomogeneity. The analogy between the two proce- 
dures also explains the formal similarity between the GP equation for condensed composite bosons at zero 
temperature and the Ginzburg-Landau equation for the superconducting order parameter near the critical 
temperature. 

Note that the boson-boson scattering length as that can be introduced in Eq. ljSTj) equals 2a p, a value 
that corresponds to the Born approximation for the scattering between composite bosons. This value is 
borne out by any mean-field treatment of the BEC limit of the BCS-BEC crossover. Improvement on the 
treatment of the low-energy scattering between composite bosons results in different values for the ratio 
o-b/clf 03]- In particular, an exact treatment [55] EH] yields the value ae/aj? = 0.6. Although it is possible 
to include the fermionic scattering processes responsible for this reduced value of as/aF in the derivation 
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of the GP equation for the composite bosons [70] (which thus would contain the correct value 0.6ap of as 
instead of 2a p), we regard this quantitative refinement in the BEC limit to be beyond the scope of the 
present treatment, where we concentrate instead in solving the BdG equations © throughout the whole 
BCS-BEC crossover. 

Note further the essential role played by self-consistency in the derivation of the GP equation, since for 
this equation self-consistency can by no means be dismissed. This is true in spite of the fact that the BdG 
equations |5J are apparently linear equations, while self-consistency is introduced by the condition ([6]) on 
the local gap parameter A(r). 

The derivation of Eq. (|5Tj) from the BdG equations provided in Ref.[SD] relies on the smallness of the 
parameter A/|/i| when 1 (kpap) -1 , as well on the smoothness condition to be satisfied by the potential 
V(r). In practice, however, even for a square- well potential with sharp edges we have favorably compared the 
results of our numerical solution of the BdG equations (J5J) for large enough values of the coupling parameter 
(kpap) -1 with the independent numerical solution of the GP equation (|51| . obtained with the same value 
of of and with the value of the bosonic chemical potential fiB = 2/x + (map) -1 expressed in terms of the 
corresponding fermionic chemical potential /x. 

An example of this comparison is shown in Fig l 101 when (kpap)^ 1 = +3.0 for rectangular and Gaussian 
barriers of comparable heights and widths. Two alternative physical situations are further considered in this 
figure, namely, when the size of the composite bosons (given by £ pa i r = ap/y/2 in this limit [12]) is smaller 
or larger than the barrier width. Only in the latter case, in fact, we expect deviations to occur between the 
solutions of the BdG and GP equations, since the composite nature of the bosons mostly manifests itself 
when they cannot fit into the barrier region. Moderate deviations along these lines are evident from FigfTUl 
We shall return to this interesting issue in Section 4 while providing an overall discussion of the results of 
the numerical calculations. 




Figure 10: Comparison of the solutions of the BdG equations (dashed lines) and GP equation (full lines) for 
(kpap)^ 1 = +3.0. The magnitude |$(a;)| is here reported for q = in units of its bulk value 4>o. Two (a) 
rectangular and (b) Gaussian barriers of widths Lkp = akp^/2~rr = (0.1,2.0) and height Vo/Ep = 4.0 are 
considered [the upper (lower) curves correspond to the narrow (wide) barriers] . The barrier widths need be 
confronted with the size £ pa ir&F = (3\/2) _1 = 0.23 of the composite bosons for this coupling. 



Besides being an important benchmark for the numerical solution of the BdG equations in the BEC 
limit, the independent solution of the GP equation plays in this case an additional role as it permits one to 
explore parameter ranges (concerning, for instance, the barrier) for which a direct numerical solution of the 
BdG equations becomes exceedingly difficult and/or time consuming. In these cases, a preliminary solution 
of the GP equation may serve as an exploratory mean for focusing on the physical effects one is after. In 
the specific context of the Josephson effect, the solution of the GP equation may also serve for evidencing 
an alternative physical meaning of the effect itself, as discussed in Appendix C. 

Still regarding the Josephson effect, the above mapping between the BdG and GP equations in the BEC 
limit allows one to establish a direct link with the literature concerning coherence effects in Bose-Einstein 
condensates, which are made up of point-like bosons and are themselves invariably described by the GP 
equation [71] . 



24 



Finally, we refer the reader to Appendix D for some useful manipulations one needs to perform on the 
GP equation (|5ip in order to ease its numerical solution in the presence of a barrier with a slab geometry 
and a finite current. 



3 Numerical strategies and procedures 

We pass now to discuss a number of procedures that prove necessary for an efficient numerical implementation 
of the method discussed in Section 2, whereby the solution of the BdG equations ([5]) in the presence of a 
barrier with a slab geometry was set up operatively. In particular, we will discuss the numerical strategies 
that are required to speed up the cycles of self-consistency for Ea. p^j) . a process that proves especially 
important when approaching the strong-coupling (BEC) limit. 

To reconcile an accurate solution of the scattering problem based on the BdG equations with a feasible 
effort for attaining self-consistency, we have considered at the ouset two different spatial grids, one grid at 
the positions {xi; I = 1, • • • , M — 1} where the continuity equations for the scattering problem are specified, 
and the other grid at the positions {xf,£ = 1, • • • , M — 1} where the local gap parameter A(x) is taken as a 
variable in the cycles of self-consistency. In practice, the grid {%} can be taken to be a coarse subset of the 
grid {xe} (such that M <C M), while the values of A(a^) between two successive points Xj> and % +1 (where 
A(x) is specified in the cycles of self-consistency) are obtained by suitable interpolation (which is necessary 
because in the scattering problem the values of A(x) need be specified over the denser grid {a;^}). Typically, 
we have taken M ~ 100 and M ~ 20. Furthermore, a judicious choice of the points {%} of the coarse 
grid proves necessary, so as to avoid obtaining sharp variations of A(x) where the barrier itself has jump 
discontinuities (this shortcoming is especially evident in the BEC limit when comparison with the results of 
the GP equation is possible). To this end, it proves sufficient not to include in this grid the points where 
these discontinuities occur. 

Equation (|29p is then solved separately in each of the seven energy ranges introduced in sub-section 2b. 
Special care is necessary in ranges 1 and 5 where bound-state solutions occur, as these may sometimes escape 
even a quite refined search. As a matter of fact, precise inclusion of these bound-state solutions turns out to 
be essential for an accurate determination of the physical quantities to be calculated, and also for reaching 
the self-consistency condition itself of the BdG equations. In addition, in range 2 resonant (quasi-bound) 
states may occur which need also to be treated with sufficient accuracy. 

Integration over k = (fcx,kii) is next performed in the self-consistency equation fD^i as well as in the 
equations ([5]) for the density and @ for the current. As discussed in Appendix A, the integration over 

need be performed first (and contains also a discrete sum over the bound states that are defined for 
each given value of k||). In this integration (which extends in principle from — oo to +oo), the variable k± 
coincides (in a piecewise continuos fashion) alternatively with one of the (real) wave vectors reported in 
Figs and O depending on the energy ranges of FigJU 

Both integrals (over k± and kn) are further split into two ranges. For each integral, a cutoff k c is 
introduced (which may, of course, differ for the _L and || directions), such that for k < k c the integration is 
done numerically while for k c < k < +oo one first extracts the appropriate power-law decay of the integrand 
together with the relative coefficient and then perform the integration analytically. Since Eqs.®, J9]), and 
(frjj are all local relations, also k c depends parametrically on x, the largest values of k c (x) occurring at the 
position (x = 0) of the barrier. 

With these considerations in mind, we plot in Fig[TT]the values of k^(x = 0) (full line) and k^(x = 0) 
(dashed line) vs the coupling parameter (kpap) -1 for a typical SsS barrier we have considered (we have found 
that these plots depend only mildly on the chosen barrier). Since fc^_(x = 0), in turn, depends on the value 
of fcy (x = 0), we have here considered k^x = 0) in the extreme case when k\\ (x = 0) = k^(x — 0). Note that 
both quantities increase rapidly when approaching the BEC limit of the crossover (where 1 <C (/cfOf) -1 ), 
thus evidencing the increasing numerical effort one has to sustain on the BEC side with respect to the BCS 
side of the crossover. 

Even though Eq. (|14[) was originally introduced as the condition for enforcing self-consistency, a "mixed 
choice" is actually to be preferred for the purpose. It turns out, in fact, that to the real and imaginary parts 
of Ea. (Ti"4")) there correspond different sensitivities to the self-consistent parameters {A(xi)} and {(p(xi)} as 
well as to the current wave vector q. Specifically, the imaginary part of Eq. (|14[) proves to be quite insensitive 
to the values of these variables, so that it cannot be the appropriate candidate for their self-consistent 
determination. The real part of Eq. (fT4"]) . on the other hand, is suitable to the purpose. We have then 



25 




-1 -0.5 0.5 1 

(k F a F ) 



Figure 11: Cutoffs k?,(x — 0) (full line) and kj_(x — 0) (dashed line) in units of kp, which are used in 
the numerical integrations, vs the coupling (kpap)^ 1 for a typical rectangular barrier with Lkp — 4.0 and 
V /E F = 0.1. 



replaced the imaginary part of Ea. (fT4")) by the current equation ([5]), with its left-hand side set equal to 
j = qno /m where no is the asymptotic (bulk) value of the density far from the barrier (the current being 
uniform). The choice of Eq.Q is also suggested by the fact that, in the strong-coupling (BEC) limit, the 
continuity equation (and thus the constancy of the current) results from taking the imaginary part of the 
GP equation [cf. Appendix D] . 
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Figure 12: Evolution toward self-consistency of |A(x)| and 2<p(x) for the coupling value (/cf^f) -1 = —0.5 
with a rectangular barrier of width Lkp = 5.0 and height Vo/Ep — 0.5. Dotted, dashed, and full lines stand 
for initial, intermediate, and final values of the functions, in the order. 



With the two above equations, the approach to self-consistency has been achieved via the multi-dimensional 
Newton method, which proves to be more efficient than the iterative method (used, for instance, in Ref. 57 ) 
in reducing the number of cycles. By the present method, self-consistency was invariably reached within at 
most five cycles. As an example, we plot in FigU^a typical evolution toward self-consistency of the profiles 
of |A(x)| and 2<j>{x). 
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The numerical calculations have been performed by imposing the uniform value qrio/m of the current, 
with the provision that this value should be self-consistently sustained. We will, in fact, find that to each 
coupling (fefOi?)" 1 there corresponds an upper (critical) value q c that can be self-consistently sustained, 
which we will associate to the Landau criterion for superfluidity as discussed in sub-section 4f. By this 
procedure, a value S<f> for the asymptotic phase difference 2[<ft(x = +00) — <p(x — —00)] between the two sides 
of the barrier results for given q. Alternatively, the numerical calculations can be performed by fixing the 
value of 4>(x = +00) — <f>(x — —00) and letting q be an independent variable to be detemined self-consistently. 
We have actually adopted either one of these two procedures depending on the circumstances. Specifically, 
when the value of q is selected in advance, we have found that the numerical solution invariably converges 
to that corresponding to the "left" branch of the Josephson characteristics (see sub-section 4b) . To access 
also the solutions corresponding to the "right" branch of the Josephson characteristics, therefore, the value 
of 8(j> need be selected in advance. The corresponding stability of the solutions of the left and right branches 
of the Josephson characteristics will be addressed in detail in sub-section 5b. 

As a matter of fact, we have found it more convenient to use the derivative <fi'(x) of the phase (and not 
the phase (f>(x) itself) taken at the coarse set of points {xj}, as the independent variables to be inserted 
in the self-consistency conditions. That this is a more appropriate choice of variables is also evident from 
the solution of the GP equation [cf . Appendix D] , although it appears less evident when solving the BdG 
equations troughout the BCS-BEC crossover since in this case (f>'(x) occurs explicitly only in the current 
equation. Operatively, we begin by an initial guess of <j/j at the points {%}, from which we obtain a 
continuous function <p'(x) by (a cubic spline) interpolation, and then obtain the required values (j>i for each 
I as follows: 

fa = / dx' cj>'(x') . (53) 
J -co 

Numerically, the advantage of this procedure over the more obvious choice of 4>(x) as the independent variable 
(from which <p'(x) could be obtained by numerical differentiation when necessary), is that the integration 
introduces a smaller numerical error than the differentiation, so that one can consider a more limited number 
(M — 1 ) of points in the coarse grid. 

An additional numerical difficulty we have encountered when solving the equation (|29|) is that some 
matrix elements (|28[) may become exceedingly large owing to the presence of the factor exp{ik^ x} therein. 
To better appreciate this point, let's refer, for instance, to an electron-like excitation impinging from the left 
in range 4. In this case one sees from the forth panel of Figf5]that k^ possesses a negative imaginary part, 
in such a way that exp{ik^x} grows without bound when x becomes large to the right of the barrier. This 
was precisely the reason for excluding the corresponding wave function (|24|) from the linear superposition 
(1231) in the outermost left interval with I = M [cf. Table 1], thus setting cm = at the ouset. The same 
assignement cannot a priori be done, however, with the coefficient cjf-i of the next interval to the right, a 
coefficient that should turn out from the solution of Eq. (j2"!)|) to be exceedingly small in order to compensate 

for the large factor expjifcj^^ir}. By this mechanism, one would inconveniently multiply quite large with 
quite small numbers, a procedure that may lead to substantial numerical errors. To reduce these errors, 
one may absorb the large factors expjifc^™ x} (whenever they occur) in the associated coefficients of the 
expressions (|25p . by redefining in each interval: 

&i(x;q,E) = a e Ty(x — xi-i;q, E) + be Ty(x — xe-i;q, E) 

+ c i rf\x~x l - l -q,E)+d l rf\x-x l - l -q,E) (54) 

where ae = ae exp{ik9~ xt-i}, be — be expfykj? X(.-i}, ce = Q exp{ikjpX£-i}, and dt = dt exp{ik^X£-i}. In 
this way, the continuity conditions (126p are replaced by 

M e {x = x t - xi-i) We = Me+i{x = 0) W e +i (55) 

where Wj — (a,f,be,ce,de) [cf. Eq. (j2"T]) ]. The expressions entering Eq. (|2"5|) are consequently modified. 
By this expedient, we have effectively enlarged the range of the energy E which is required for accurate 
numerical calculations (especially in the BEC limit of the crossover). This is because the phase factors 
exp{ik^ x} entering the matrix elements (|2"8)l can blow up because either Imjfc^} or |x| get large. With 
the transformation ([5"4")l . \x\ becomes at most limited by the width of the largest interval Axe = xe — xe-i- 
Typically, Ax^Im{fc} max <~ 500 in order to avoid numerical under- or over-floating. Here, one can choose 
Axe ~ £ P hasc/10 in order to obtain the profile of the gap with sufficient accuracy (apart from the occurrence 
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of the Friedel oscillations over a scale kp in the BCS limit of the crossover, where smaller values of Im{fc} max 
need be considered). This implies that £ p hasc Im{fc} max ~ 5000. Since (kp £ p hasc)max ~ 10 in the relevant 
range of the crossover (cf. Fig[7|), one obtains lm{k} max /kp ~ 500 as a conservative estimate for the upper 
bound on Im{fc} max while Axg kp ~ 1, values that turn out to be totally sufficient for all practical purposes. 

Note further that the present method does not put too severe a limitation on the maximun number of 
zones M one has to rely on for an accurate description of the problem (the major limitation originating, in 
practice, by the increasing computational time when M increases). This contrasts with the transfer- matrix 
method described in Appendix E, where the number of zones is severely limited by the numerical under- or 
over-floating problems one encounters when M exceeds 40. 



4 Results for an SsS barrier 

We pass now to give a detailed account of the results that we have obtained by solving numerically the 
BdG equations © subject to the self-consistent condition © in the presence of a stationary superfluid flow 
impinging on a barrier, according to the methods described in Sections 2 and 3 which enable us to scan 
the whole BCS-BEC crossover at zero temperature [IS]. Our systematic calculations will not only provide 
a bridge between the treatments of the stationary Josephson effect which were so far considered separately 
for fermionic and bosonic systems in the literature, but they will also uncover several novel and interesting 
features about the Josephson and related effects that could not emerge when addressing these (BCS and 
BEC) limits separately. 

4a. Spatial profiles of the fermionic wave functions and the complex gap 

parameter 

Knowledge of the solutions u„(r) and tv(r) of the BdG equations © is required to calculate the gap 
function ^ as well as the number density ([5]) and current © (or the total energy - cf. Eq. (f8T)]) of sub-section 
5a). Even though these wave functions have no individual meaning, it is nevertheless instructive to examine 
their spatial behavior along r = (x, 0, 0) in some characteristic cases to get a physical feeling about their 
overall role. 

As an example, in Fig[T3] we plot |wj/(a:)| 2 and |iv(a;)| 2 vs x at self-consistency for the following cases 
that correspond to two of the four ranges 




xk F xk F 

Figure 13: Plot of |it„(:z;)| 2 (full lines) and |u„(x)| 2 (dashed lines) vs x at unitarity for (a) a bound-state 
solution and (b) an electron-like excitation in the continuum, for a rectangular barrier in the presence of a 
finite current. The values of parameters are specified in the text. 



identified in Fig|4ja) : (a) A bound-state solution occurring in range 1; (b) An electron- like excitation whose 
energy E is comprised in range 3. These states are calculated for a rectangular barrier of width Lkp = 4.0 and 
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height Vo/Ep — 0.2, with k|| = and q/kp — 0.128, and for the coupling value (Ufclf) 1 = at unitarity. 
Both the continuum and the bound-state solutions are affected by the depression of the magnitude |A(x)| of 
the gap parameter about x = 0, which originates from the presence of the repulsive barrier. [In particular, 
the bound state would not exist were not for the depression occurring in |A(x)|.] 
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Figure 14: The profiles of |A(x)| and 2<j){x) vs x are shown for the same parameters of Fig ll3l The arrow 
identifies the extension of the length £ p hase- 



The corresponding profiles of the magnitude |A(x)| and phase 2<f>(x) are reported in Fig[l4]for the same 
parameter values. [Recall that, by our definition Eq. fTS)) . the phase of the local gap parameter is 2cf>(x) and 
not <j)(x).] Note how the spatial variations of these quantities occur over the scale £ p hase introduced for the 
homogeneous case (as identified by the arrow in the figure), and that these variations are mutually related 
by the condition of the current being uniform. 
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Figure 15: Evolution of the profiles of |A(x)| and 2<f>(x) from weak to strong coupling, obtained for a 
rectangular barrier with Lkp — 4.0 and Vo/Ep = 0.1 (the values of the coupling are specified in the text). 
The inset shows the corresponding density profiles n(x). The solution of the GP equation for the coupling 
value on the BEC side is shown for comparison (full line). 



We will devote sub-section 4e below to a detailed analysis of the bound-state solutions (the so-called 
Andrccv-Saint James states), since they play a special role in the theory of the Josephson effect. We 
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examine instead here the way the profiles of |A(a;)| and 2<p(x) for given barrier and q evolve with increasing 
coupling, by considering in Fig[l5]three characteristic coupling values (kpap)^ 1 — —0.8 (dashed lines), 0.0 
(dotted lines), and +3.0 (dashed-dotted lines), for which fcF^phasc = (1.33,0.68, 1.32) from weak to strong 
coupling, in the order. In each case, the value of q corresponds to the maximum sustainable value of the 
current for the given barrier and coupling. The profiles obtained by the independent solution of the GP 
equation (|5ip are also reported for comparison (full lines) when (kpap) 1 = +3.0, showing a remarkable 
agreement between the two independent (BdG and GP) calculations for this case. 

Note in Fig[15] the progressive depression of |A(x)| at the center of the barrier as the coupling is in- 
creased, with an associated steepening of 2<p(x) that results in an increase of the asymptotic phase difference 
5(f) = 2[4>{x = +oo) — (j>(x = — oo)] between the two sides of the barrier. The depression of |A(x)| is also 
reflected in the corresponding density profiles shown in the inset (where no is the bulk fermionic density). In 
particular, |A(x)| 2 and n(x) are proportional to each other in the BEC limit (cf. Ea.Q52])). These features 
are independent from the shape of the barrier, as we have explicitly verified by repeating the calculation 
using a Gaussian barrier with the same area for definiteness (that is, taking a = L/^/2%). 



4b. Josephson characteristics: Current versus phase relations 

The above calculations, which relate the value of the asymptotic phase difference 5<j> to the value q of the 
current wave vector, can be performed in a systematic way so as to obtain the function q(S(f>) which identifies 
the characteristic Josephson current-phase relation J(6(f>) where J = qn^jm. In this context, one actually 
finds that two different values of 6(j> (corresponding to two different spatial profiles <fi(x)) are associated with 
a given value of q. These two solutions give rise to two different branches of the J vs <p relation, which 
we shall conventionally call the "left" and "right" branch and which will be shown in Section 5 to possess 
different stability properties against the introduction of fluctuations. Setting apart for the moment the issue 
of their stability, we will represent here the two branches as a single continuous curve and discuss their 
overall properties across the BCS-BEC crossover. 

In Fig|16l we report the behavior of J vs 6<j) for given barrier and different coupling values spanning 
the crossover region. A rectangular [Fig lTCT a)] and a Gaussian [Fig lTBT b)] barrier are considered with the 
same value of the (dimensionless) area Zp = mZ/kF (which equals {Lkp){yo / Ep) /2 for a rectangular and 
(akF)(Vo/EF)y/w/2 for a Gaussian barrier - cf. Ea.([39])). 
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Figure 16: Josephson characteristics q(54>)/kF = J{Scf))m/ (fcjmo) f° r ( a ) a rectangular barrier with Lkp = 4.0 
and (b) a Gaussian barrier with (akp) = 1.6 (Vo/Ep = 0.1 in both cases). The three curves correspond to 
(kFdF)^ 1 — —1.0 (dotted lines), 0.0 (dashed lines), and +1.5 (full lines). 



Two features appear evident from these plots. In both cases, the Josephson current is considerably 
enhanced at unitarity with respect to the BCS and BEC sides (a phenomenon to which we will return more 
extensively in the following). In addition, the curves are seen to stretch from being proportional to sm(5<fi) 
to being (almost) proportional to cos(<5</>/2) following the evolution from strong to weak coupling. 
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The proportionality to sin(50) when approaching the BEC side with given barrier height represents a 
standard result for the Josephson characteristics, and reflects the fact that in this limit the barrier height 
Vo becomes larger than the chemical potential [Ib of the composite bosons, which is the only other relevant 
energy scale of the problem. 

The proportionality to cos((50/2) when approaching the BCS side, on the other hand, results from having 
achieved full self-consistency in the calculation, whereas the non-self-consistent calculation discussed in sub- 
section 2c (cf. Ea. (|50|) ) would yield a sin(50/2) proportionality for a weak barrier with Zp <C 1 (Zp = 0.2 
for the barriers of Fig [16]). 
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Figure 17: Josephson characteristics q(5<j))/kp = J(6(f))m/(kpno) for (a) a rectangular and (b) a Gaussian 
barrier with the same widths of Fig [11)1 but with the larger height Vo/Ep = 0.4. The coupling values are the 
same of FigJTHl 



In FigfTrjI a rather small value was considered for the ratio Vo/Ep. [Note, however, that such a value for 
the height of the barrier becomes soon comparable with the chemical potential \ib of the composite bosons 
when approaching the BEC side of the crossover, and in this sense it is not small.] In Fig[l7|we report for 
comparison the corresponding behavior of J vs 8(j> for the larger value Vo/Ep = 0.4, the other parameters 
remaining the same as in Fig |16l In this case, all curves are much suppressed with respect to those of Fig |16l 
and approach the sin(<5</>) dependence that characterizes a strong barrier not only on the BEC but also on 
the BCS side of the crossover. 

Conversely, in Fig ll8l we show how the Josephson characteristics evolve when varying the strength of a 
rectangular barrier while keeping the coupling fixed at unitarity. The example shown is for a rectangular 
barrier with height Vo/Ep spanning a wide range. One observes again the progressive evolution, from 
a sm(5<fi) dependence for strong barriers to an approximate cos((5</>/2) dependence for weak barriers [72] , 
We have found that this behavior is not restricted to unitarity, but is rather generally reproduced for any 
coupling. 

Before discussing the additional features of the Josephson characteristics over the whole BCS-BEC 
crossover, we restrict ourselves for a moment to the peculiar effects that can be identified from the nu- 
merical calculations in the two extreme BCS and BEC limits. 



4c. The weak-coupling (BCS) limit and the occurrence of Priedel oscillations 

We have already mentioned that for given coupling the value of £ p hasc sets the scale of the upper limit 
for the spatial range that need be explored by our numerical calculations about the position of the barrier, 
because only past this distance the magnitude |A(x)| of the gap parameter recovers its bulk value Ao. In 
Section 3, we have accordingly estimated the value 10 as a resonable upper bound for the product fcF^phase, 
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Figure 18: Evolution of the Josephson characteristics q{8<f>)/kp = J (64>)m / '(hp no) at unitarity for varying 
barrier height. Three different values of Vo/Ep = 0.025 (clotted line), 0.1 (dashed line), and 0.4 (full line), 
are considered for a rectangular barrier of width Lkp = 4.0. 

which in turn implies from Fig[7]that the coupling value (kpap)^ 1 — —2 represents, in practice, the limit as 
far as one can progress toward weak coupling. 

In Fig ll9l we report a typical profile of |A(x)| in weak coupling, obtained for (kpap)^ 1 = —1.5 when 
^F^phasc = 3.4. In this case, a supercurrent with q/kp = 0.01 impinges on a rectangular barrier of width 
-£/£phase = 0.01 and height Vo/Ep = 3/(Lkp) that resembles a delta function. 

Note from this figure the presence of Friedel oscillations with characteristic wave vector 2kp, which 
originate from the sharpness of the Fermi surface in wave- vector space and distintictly modulate the profile 
of |A(x)|. Apart from these oscillations, note also the presence of two length scales in the overall profile of 
|A(x)|, as discussed in sub-section 2c: A first scale (given by kp l in this case) within which |A(x)| is strongly 
suppressed locally by the barrier, and the second scale £ p hase over which |A(x)| smoothly approaches its bulk 
value A . [The additional lenght scale Cy — (2mVb)~ 1 ^ 2 associated with the barrier does not affect the 
profile of |A(x)| in the present case, since kp£v — 0.106.] 
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Figure 19: Profiles of |A(x)| when (kpap) 1 = —1.5 and qjkp = 0.01, for a rectangular barrier of width 
Lkp = 3.4 x 10~ 2 and height V /E F = 88.24. 
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In weak coupling, we have thus far considered barriers resembling a delta function for which L/£ p hasc 1. 
It is, however, interesting to consider also barriers whose width L is comparable with (or possibly even larger 
than) £ ph 
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Figure 20: (a) Profiles of |A(x)| for (kpap)^ 1 = —1.5 and q/kp = 0.01, obtained with three rectangular 
barriers having Lkp = 0.01 (full line), Lkp = 2.6 (dashed line), and Lkp = 5.3 (dotted line), while Z — 0.525 
in all cases, (b) Comparison of the profiles of |A(x)| (full line) and n(x) (dashed line) for the narrowest 
barrier considered in (a). 



The profiles of |A(x)| obtained for the coupling (kpap)^ 1 = —1.5 and current q/kp = 0.01, and cor- 
responding to three different values of L which are respectively smaller than, comparable with, and larger 
than £ p hase, are shown in Fig l20f a). This figure makes it evident that there arc definite differences between 
narrow and wide barriers, and that the length scale £ p hase is responsible for these differences. 

The profile of |A(x)| for the narrowest barrier of Fig |20l a) is further compared in Fig |20l b) with the 
number density n(x), showing a characteristic de-phasing of the Friedel oscillations affecting these two 
quantities. This implies that |A(x)| and n(x) are somewhat correlated with each other even in the weak- 
coupling limit (while in strong coupling we know from the relation (|52[) that |A(x)| 2 is strictly proportional 
to the bosonic density). 

4d. The strong-coupling (BEC) limit and the internal structure of the 

composite bosons 

We have already discussed at some length in sub-section 2d the role played by the bosonic GP equation 
when the fermionic BdG equations are solved in the BEC limit of the BCS-BEC crossover, that is, when 
fermion pairs of opposite spin bind together so as to form composite bosons whose size is much smaller than 
the average interparticle distance. 

Even in this limiting case, however, the presence of a barrier may affect the comparison between the 
results of the corresponding solutions of the GP and BdG equations, when the size of the composite bosons 
is comparable with the width of the barrier. Under these circumstances, the composite nature of the bosons 
plays a role and the fermionic nature of their constituents emerges, when a detailed comparison is made with 
the solution of the GP equation which completely ignores the internal structure of the bosons. 

To spot these kinds of effects with more evidence than in Fig[TUl we report in Fig|2T] the self-consistent 
profiles |A(x)| for the lesser extreme coupling value (kpap) _1 = +1.5 and for a rectangular barrier with 
two different widths Lkp = (0.1,4.0) and heights Vo/Ep = (0.1,0.4). For this coupling value, the size of 
the composite bosons equals /cF^pair = 0.44, and it is respectively larger (upper panels) and smaller (lower 
panels) than the barrier width for the cases here considered. In each case, the solution of the GP equation 
is also shown for comparison. 
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Figure 21: The profiles |A(x)| when (^pcti?) -1 = +1.5 (dashed lines) are compared with the corresponding 
solutions of the GP equation for the same coupling (full lines), for four rectangular barriers with Lkp = 
(0.1,4.0) and Vq/Ef = (0.1,0.4). In all cases here considered, we have fixed the value Sip/ir = 0.1. 

The comparison evidences definite deviations between the two (BdG and GP) calculations in the region 
close to the barrier when the size of the composite bosons is larger than the barrier width, such that the 
internal wave function of the fermion pair is considerably affected by the presence of the barrier. Specifically, 
taking into account the composite nature of the bosons as in the BdG equations leads to a noticeable 
depression |A(x)| in the vicinity of the barrier with respect to point-like bosons. The panels of Fig|2T1make 
it further evident that the height of the barrier has little effect on the degree of discrepancy between the 
BdG and GP calculations. 



4e. The role of the Andreev-Saint-James bound states 

The Andreev-Saint- James reflection was originally conceived as occurring at an interface between a 
superconductor and a normal material [51] . By this process, an electron incident on the interface from 
the normal material (with an energy smaller than the superconducting gap) is retro-reflected as a hole of 
opposite spin and momentum, while a Cooper pair forms on the superconductor side. This occurs because 
an electron impinging from the normal side at energies smaller than superconducting gap cannot propagate 
inside the superconductor, while Cooper pairs can do so. At the same time, the reflected hole ensures current 
conservation. 

This process is actually a key feature for the solutions of the BdG equations. As such, it also occurs in 
the different situation considered in this paper when a barrier separates two superconductors. In this case, 
the multiple Andrccv-Saint- James reflections occurring at the barrier sides may interfere constructively to 
give rise to Andreev bound states localized at the barrier. These bound states, which emerge when solving 
the BdG eigenvalue problem even with a repulsive potential barrier, originate from the lowering of the gap 
profile at the potential barrier, which acts as an effective binding potential in the off-diagonal terms of the 
BdG equations. 
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The importance of these bound states with subgap energies stems from the fact that, following the 
pioneering work of Kulik [73] ■ they are believed to provide the main channel (pictorially, a bridge) for 
Cooper-pair transfer via Andreev-Saint-James scattering between the two superconductors separated by a 
barrier 64, 52J. As a consequence, the process of Andreev-Saint-James scattering is intimately related to the 
dc Josephson effect when a pair of correlated fermions is transferred from one to another superconductor [3] . 
It is also worth mentioning that the spectroscopy of these levels constitutes a quite active field of research 
both in conventional as well as in cuprate superconductors [63] . 

The above considerations, about the connection between the Andreev-Saint-James reflections and the 
dc Josephson effect, have historically been applied to weak-coupling situations where the presence of the 
Cooper pairs play a major role. With the present calculation based on the self-consistent solution of the 
BdG equations, whereby the fermionic attraction can be varied at will from weak to strong coupling, we are 
in a position to address the connection between the Andreev bound states and the flow of the Jospehson 
current within a more general context, up to the point when the Cooper pairs give the way to the presence 
of composite bosons whose underlying fermionic degrees of freedom become progressively more irrelevant 
for increasing fermionic attraction. In addition, our need for taking into account the fermionic wave vector 
transverse to the current flow (in order to be able to form the composite bosons in the strong-coupling limit) 
brings out the issue of the dispersion of the Andreev bound states as the magnitude of this wave vector is 
varied. 
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Figure 22: Fraction of the superfluid current carried by the Andreev bound states vs (A^of) -1 at the position 
of a barrier with Lk F = 4.0 and V /E F = (0.05, 0.10, 0.20, 0.40) from bottom to top. 



We have already discussed in sub-section 2c that, in the weak-coupling limit, a non-self-consistent calcu- 
lation suffices to obtain a reliable description of the Josephson effect for a Dirac-delta barrier whose strength 
is large enough, and that in this case a single Andreev bound state accounts for the whole Josephson current 
at the position of the barrier (a result that was well-known in the literature - cf. Ref.[74]). The point is here 
to understand to what an extent the above statement survives to (i) the introduction of self-consistency, (ii) 
the presence of a finite barrier width, and (iii) the evolution toward strong coupling. 

To this end, wc plot in Fig |22l the ratio of the superfluid. current jAndreev 

carried at the position of the 

barrier through the Andreev bound states with respect to the total current j, across the BCS-BEC crossover. 
The results are shown for several barriers and include the contribution of the whole Andreev band spanned 
by the wave vector k|| transverse to the current flow. Note how the contribution of the Andreev bound states 
may sometimes exceed 100% of the whole current; in this case, the continuum levels contribute a current 
flowing in the opposite direction (this feature was already pointed out in Ref. 54J for the weak-coupling 
limit). 
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Figure 23: Square magnitude of the Andreev wave function with k|| = vs xkp across the position of a 
barrier with Lkp = 4.0 and Vq/Ef — 0.10, for three coupling values: (fc^a^) -1 = —1.0 (full line), 0.0 
(dashed line), +1.0 (dotted line). In each case, the superfluid current is taken at its maximum sustainable 
value. 



In addition, Fig|23l shows the degree of localization of an Andreev bound state with ky = 0, by plotting 
the square magnitude |u^(a;) | 2 + |u„(a;)| 2 of the corresponding wave function vs x for three characteristic 
coupling values across the BCS-BEC crossover. One sees from this figure that the spatial extension of the 
Andreev bound state decreases for increasing coupling from the BCS to the BEC limits. 

We shall return in sub-section 4g to the special role played by the Andreev bound states on the BCS side 
of the crossover with respect to the stability of the superfluid flow. 

4f. Critical Josephson current and its relation with the Landau criterion for 

superfluidity 

We have seen in sub-section 4b that the overall shape of the Josephson characteristics varies considerably 
by changing the width L and height Vq of the barrier as well as the value of the coupling throughout the 
BCS-BEC crossover. In particular, these variations affect the maximum Josephson current, which we have 
seen to increase as the ratio Vq/Ef is decreased. Since the maximum Josephson current is probably the most 
important feature of the Josephson characteristics for practical purpose, it is relevant to address the question 
whether there exists an intrinsic maximum allowed value of the Josephson current for any given coupling. In 
other words, we shall address the question whether adopting the limiting procedure of progressively lowering 
the barrier height leads to an intrinsic upper value of q c /m for the maximum velocity of the superfluid flow 
at any given coupling. 

One may anticipate on physical grounds that a vanishingly small barrier acts as an impurity that probes 
the stability of the homogeneous superfluid flow. In this sense, it plays a similar role to the walls of the 
container in the context of the Landau criterion for superfluidity [53] . One thus expects that the velocity 
q/m of the superfluid flow never exceeds the corresponding value of the critical velocity as obtained by the 
Landau criterion. This critical velocity is determined by the onset of quasiparticle excitations which are of 
a different nature on the two sides of the crossover, namely, pair-breaking excitations on the BCS side and 
sound-mode quanta on the BEC side. 

Before showing in detail how our numerical self-consistent solutions of the BdG equations indeed capture 
this important physical feature on both sides of the crossover, we shall briefly consider the mechanism 
through which the BdG equations are able to recognize the occurrence of an instability when the velocity of 
the superfluid flow increases beyond a certain limit. 

We begin with a preliminary discussion about the role of Galilean invariance in the BdG equations. We 
thus consider how the wave functions Uk(r) and «k(r), that are solutions of Eq.([5]) in the absence of a barrier 
(homogeneous case), are transformed when passing from the lab frame if to a frame K' which moves with 



3G 



velocity V with respect to K. This transformation reads: 



Wk(r) \ _ iq . r / Uk'(r) 
«k(r)* J ~ I «k'(r)* 



(56) 



where q = mV, k = k' + q, Uk'(r) = Uk' e lk r , and i>k'(r) = i>k' e lk r (with ttk' and v-^i real). Solving 
for the BdG equations alternatively in the two frames then yields the following relations between the gap 
parameter, chemical potential, and energy eigenvalues in the two frames: 

A(r) = e 2iqr A' , £* = // + ^— , e k = «£, + ^ • k' (57) 

zm to 



where e k/ = y (k' /(2m) — /i') 2 + A' 2 . These relations have been already utilized while discussing Eas. lfTS")) - 
(fT7|) near the end of sub-section 2a, although with a slightly different notation. Note that in the above 
relations, no factor appears that could limit the value of the relative velocity V. For a truly homogeneous 
superfluid, this feature would naivly reflect itself in the absence of an upper value for the current j = nV 
flowing in the frame K at rest with respect to the superfluid, in contrast to the Landau criterion. 

To appreciate the way the BdG equations recognize the occurrence of an instability when the velocity of 
the superfluid flow increases beyond a certain limit, we rewrite the gap equation ^ for the homogeneous 
case in the form: 

1 r dk' 1 

g ~ J (2^^ 

with reference to the newly introduced notation in the frame K' where the superfluid is at rest. Consider, 
in particular, an infinitesimally small (albeit non-zero) temperature for which the presence of the Fermi 
function can be neglected whenever its argument is positive, resulting in a gap parameter A' that does not 
depend on q. A different situation arises, however, as soon as e k , -I- — • k' vanishes, a condition which leads 
to an instability of the system felt by the BdG equations. To study this situation, we take q = (q, 0, 0) and 
set ft = fi' — [(k' y ) 2 + (k' z ) 2 )]/(2m), so that the above condition becomes: 



1 - 2/ F (4, + -i-k' 

TO 



(58) 




(59) 

where k! = k' x . Two real solutions exist for this equation provided jj' 2 + A' 2 < (p! + q 2 /to) 2 . The minimum 
value q c of q for this to occur is given by: 

^ = V^' 2 + A' 2 - n' . (60) 

TO 

This result is consistent with the Landau criterion for the collapse of the superfluid flow in the case when 
the relevant quasiparticles consist of pair-breaking excitations. 

It is also interesting to monitor the collapse of the superfluid flow directly in terms of the absorption 
of energy by the system, which becomes possible when quasiparticles can be excited. To this end, instead 
of considering the superfluid moving with velocity V with respect to a fixed barrier, we pass to the frame 
where the superfluid is at rest and the barrier moves with velocity —V. In this way, the barrier acts on the 
superfluid as an external time- dependent potential of the form 

V ext (r,t) = V(r + \t)F(t;T ) . (61) 

Here, V(r) is the static potential entering the BdG equations © while the function F(t; Tq) accounts for the 
gradual switching on and off of the barrier within the time interval (— Tq,+Tq), where To may depend on 
the velocity V. We could, for instance, take F(t;To) of the Gaussian form exp{— t 2 /^}, such that at time 
t = the full barrier appears centered at r = in an otherwise homogeneous superfluid. Correspondingly, 
the time-dependent Hamiltonian 

H'(t) = jdvp{v)V cxt {v,t) (62) 
acts on the system through the total-density operator p(r). 
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We are interested in the rate at which the external agent (represented by the moving barrier) does work 
on the system. From standard expressions of linear-response theory [75] , the total work done by the moving 
barrier is given by: 



AW = 



dp 



+ 00 



du> 
2^ 



V ext (p, u)* w Xppiv, V ext (p, w) 



(63) 



Here, x pp (p,cj) is the imaginary part of the Fourier transform of the density-density correlation function 
and K>xt(p, u>) the Fourier transform of the external potential ([GT]) : 



Wp,w) = V(p) F{lu + p • V; T ) 



(64) 



With F(t;T ) of a Gaussian form, F(oj;T ) = |T o |0r exp{-w 2 T 2 /4}. In particular, in the limit T — > +00 
one has F(w + p • V; To) — > 27r(5(w + p ■ V) in Eq. ([6"4"]) for the range of wave vectors p represented in V(p). 
This implies that the total work (|6"3"|) vanishes unless X PP {Pi w = P ' V) is different from zero for (at least) 
one of those wave vectors. 

The simplest form for Xp P (P;w) that can be studied in the broken-symmetry (superfluid) phase corre- 
sponds to the BCS approximation, which reads in the zero-temperature limit: 



X PP (P,u 



,BCS 



dk 



(u(k + pMk) + v(k + p)u(k)y 

6(E(k + p)+ E(k) +uj) - S(E{k + p)+ E(k) - u) 



(2tt) 3 



(65) 



where u(k) and u(k) are the coherence factors and E(k) the eigenvalues of BCS theory in the frame where 
the superfluid is at rest [76] (which were written as E^ in Eqs.{T]) and ([2])). Within this approximation, 
therefore, the superfluid at zero temperature is perturbed by the action of the moving barrier provided 
E(k + p) + E(k) = ±u> = =Fp • V for the range of wave vectors p represented in V(p). This implies that 



E(p - k) ± (p - k) ■ V + E{k) ± k ■ V = , 



(66) 



which is equivalent to the Landau criterion for BCS quasiparticle excitations, and requires that p = ±2k for 
the wave vector k at which a pair of these excitations begins to appear. 

Additional types of excitations come into play when considering more refined forms of the density-density 
correlation function. In particular, the BCS-RPA approximation adds the following RPA contribution [77] 
to the BCS approximation considered above for the density-density correlation function: 



Xpp(p) 



RPA 



(C(p),C(-p)) 



-A(-p) 
B{p) 



B{p) 
-A(p) 



A(p)A(-p) - B{pY 



C(p) 
C(-P) 



(67) 



where p now stands for the four-vector (p, oj). In the zero-temperature limit, the quantities A(p), B(p), and 
C{p) in Eq. (foT|) are given by the expressions [4"6] : 



A( P ) = 



dk 



(2tt) 3 



t (k + p) 2 w(k) 2 



ir\ - E{k + p) - E(k) 



v(k + p) 2 v{kf 



uj + ir] + E(k + p) + E(k) 
where 77 is a positive infinitesimal and the regularization ^ has been used; 



in 
k 2 



B(p) 



dk 



1 



i V + E{k + p) + E(k) 



ir, - E(k + p) - E(k) 

u(k + p)v(k + p)u(k)v{k) 



(68) 



(69) 



and 



C(P) = 



dk 



w(k + p) 2 



oj + ir} - E(k + p) - E(k) 
v(k + p) 2 



(70) 



iri + E{k + p) + E(k) 



u(k)v(k) 
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Note that, when taken individually, the imaginary parts of the expressions (155]) . and ([70)) arc 

nonvanishing whenever E(k + p) + E(k) = ±u>, thus leading again to the condition (|66p . When considering 
the imaginary part of the expression (|67|) . on the other hand, a novel contribution originates from the zeros 
of the denominator therein, which results in sound-mode quanta being the relevant quasiparticle excitations 
with dispersion u> = s\p\ where s is the sound velocity. In this case, the Landau criterion for the collapse of 
the superfluid yieds the well-known condition |53| 

- = a (71) 
m 

in the place of Eq.(|00|). The sound velocity s resulting from the zeros of the denominator in Eq. (p7|) 
corresponds to the Bogoliubov- Anderson mode obtained within the BCS-RPA approximation, which can be 
calculated throughout the BCS-BEC crossover [41)] . 

The two curves showing the critical wave vector gLandau vg ^ Fap ,)- 1 ) as obtained independently from 
the conditions dSUJ and ((TTJ), are plotted in FiglMlfrom the BCS to the BEC limits. In this fi gure, the "left" 
branch corresponds to the expression (|60|) appropriate to the pair-breaking excitations characteristic of BCS 
theory, while the "right" branch corresponds to the expression (|7ip appropriate to the Bogoliubov- Anderson 
mode for sound-mode excitations. In particular, in the BEC limit the value of q c of the "right" branch 
can be related to the chemical potential [Ib for composite bosons introduced in Eo. (|5ip via the relation 
q^jm — in agreement with the Bogoliubov theory for point-like bosons. 




Figure 24: The critical wave vector g^ andau at which the superfluid flow becomes unstable is shown as a 
function of (kpap)^ 1 throughout the BCS-BEC crossover. The "left" branch (blue curve) is associated with 
the appearance of pair-breaking excitations according to Eq. (|60[) . while the "right" branch (red curve) is 
associated with the appearance of sound- mode excitations according to Eq. ([7T]) . The shaded (green) area 
evidences the region of allowed superfluid flow, where the superfluid wave vector q lies below both curves. 



Note in Fig |24l that the two Landau branches cross each other near the unitary limit when (kpap) 1 = 0, 
where the maximum value of q c /m is then achieved. That superfluidity is most robust against suppressing 
mechanisms in the crossover region was pointed out in Ref.[7H] (see also Refs.[70] and [50]). 

The two branches drawn in FigJM] provide an upper limit for the the velocity V of the superfluid flow, 
because its magnitude |V| can never exceed the smallest value of the ratio q c /m corresponding to the two 
curves. The motion of the superfluid can thus occur, in principle, only within the shaded area of this figure. 

The upper value for q c /m is intrinsic in nature, in the sense that it results through a limiting procedure 
of a vanishing barrier. This is because an infinitesimal barrier acts as a "seed" for the occurrence of the 
superfluid collapse via nonlinear effects, while for a truly homogeneous system no upper limit for |V| would 
exist by Galilean invariance, as we have already emphasized. 
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The question that arises specifically when solving the BdG equations © up to self-consistency, is whether 
the curves for q c vs (kpap)~ 1 , which are obtained numerically from the maximum value q c of the Joscphson 
characteristics in the presence of a finite barrier, tend from below to the intrinsic boundary drawn in Fig |24l 
when the barrier height is progressively decreased for any value of (fcpai?)" 1 throughout the BCS-BEC 
crossover. 

That this could be a rather subtle point results by considering that, from the above procedure, one would 
naively expect to recover at most the "left" branch of Fig[24] corresponding to pair-breaking excitations, by 
arguing that the BdG equations correspond to the BCS mean field which explicitly contains only this type 
of excitations. As we have discussed in sub-section 2d, however, the BdG equations are able to recover the 
GP equation in the BEC limit. As a consequence, on the BEC side also their excitations should correctly 
correspond to the sound-mode excitations of the GP equation, since these excitations result within linear- 
response theory in terms of the BCS-RPA approximation ([57]) (cf. Ref.[77]), as confirmed by the present 
calculation. On physical grouds, this quite remarkable result has to be expected from the fact that, in the 
presence of a non-trivial geometry, the imprint of the excitation spectrum is found already in the ground- 
state wave function [45] , which in the present context corresponds to the solutions of the BdG equations at 
zero temperature for given supercurrent. In addition, reaching the upper value q c for a given barrier and 
coupling will be related in Section 5 to the onset of the energy instability of the self-consistent solutions of 
the BdG equations against the inclusion of fluctuations. 
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Figure 25: The maximum velocity q c /m obtained from the maximum of the Josephson characteristics 
is shown as a function of (kpap)^ 1 throughout the BCS-BEC crossover for several rectangular barri- 
ers, with: (a) Lkp = 2.65 and V /E F = (0.02,0.10,0.40) from top to bottom; (b) Lk F = 5.30 and 
Vo/Ep = (0.01,0.05,0.20,0.50) from top to bottom. [The width Lk F = 2.65 has been conventionally 
chosen to correspond to the value of 2kp£, P hase when (kpap)^ 1 = +3.0.] The two limiting branches of Fig|2"4l 
for vanishing barrier are also reported for comparison. 



Figure [22] shows the maximum velocity q c /m vs (kpap)^ 1 throughout the BCS-BEC crossover, as ob- 
tained from the maximum of the Josephson characteristics which result by solving the BdG equations ([5]) 
up to self-consistency for a given barrier. Values of qc/m corresponding to several barriers are reported in 
the figure. At a given coupling, these values are seen to approach from below the corresponding maximum 
values of Fig |24l related to the Landau criterion as the barrier height is progressively decreased. 

Two additional features are evident from Fig |25l First, the maximum of each curve corresponding to a 
given barrier height shifts from about unitarity toward the BCS side upon increasing the barrier height, thus 
implying that bosonic features get amplified by stronger barriers (this point will be discussed further in the 
next sub-section). Second, the rate of approach to the limiting Landau curves differs on the two sides of the 
crossover, being apparently faster on the BCS side with respect to the BEC side. 



To quantify these features, we report in Fig l26l the exponent Q of the relation 

/" Landau \ / t/ \ C 

\Qc —1c) _ I Vq ^ 



^Landau 



(72) 
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Figure 26: Exponent £ of the relation ([72]) vs (kpap) 



which is determined over the limited range 0.005 ~ Vq/Ef ~ 0.05, where T is a prefactor (this range of 
small values of Vq/Ef can also be explored experimentally, as done in Ref . [40] ) . The exponent ( turns out to 
change from the limiting value 1 on the BCS side when (kpap)^ 1 ~ — 1, to the limiting value 1/2 on the BEC 
side when 1 ~ (kpap)~ 1 , through a rather quick variation across the unitary limit (kpap) -1 = 0. These 
results imply that the present calculation, by extracting the values of q c from the self-consistent solutions 
of the BdG equations, is not only capable of recovering from below the limiting Landau curves reported in 
Figl25l but is also able to calculate the exponent that governs the approach to these limiting curves when 
the height of the barrier becomes progressively smaller. 

It is relevant to mention in this context that, after publication of Ref.[49j where the theoretical predictions 
about the critical velocities throughout the BCS-BEC crossover (like those shown in FigfJU) were reported 
for the first time, critical velocities of an ultracold supcrfluid Fermi gas throughout the BEC-BCS crossover 
have indeed been observed [40], albeit by the use of a different geometry. These critical velocities were 
experimentally determined from an abrupt onset of dissipation when the velocity of a moving one-dimensional 
optical lattice was varied. In spite of the different geometry, a pronounced peak of the critical velocity was 
indeed found at unitarity (cf. Fig. 2 of Ref.[2D]), thus confirming the physical prediction made in Ref. [45] 
about superfluidity being most robust at unitarity. In addition, the experimental shape of the critical velocity 
curve vs. (kpCLF)~ l shows the behavior predicted theoretically, of being sharper on the BCS side than on 
the BEC side of the crossover. 

In Fig. 4 of Ref. [10] the critical velocities at unitarity were further reported for varying lattice depth up 
to the point that Vq -C Ep, a limiting situation for which differences in the specific geometry should become 
irrelevant. From that figure, in the limit Vq <C Ep one can extract for q c /kp the value 0.25 and for the 
exponent £ of Eg. ([72"]) the value 0.66 (plus the value 1.23 for the prefactor T). The experimental value 0.25 
for q c lkp is somewhat smaller than the theoretical prediction 0.397 obtained at unitarity for g£ andau /fcF (cf. 
Figl2"4"]). The reason for this difference is due to the fact that this theoretical value has been obtained from 
the condition ([SU]) using the BCS mean-field values for the gap parameter and chemical potential at zero 
temperature (as reported in FigfT]). A more refined theoretical treatment would require the inclusions of 
fluctuations beyond mean field to calculate these quantities. For instance, the inclusion of pairing fluctuations 
[44] yields at unitarity the values 0.455 for fi/Ep and 0.530 for A/Ep at zero temperature (in the absence of 
a superfluid current), to be contrasted with the values 0.590 and 0.686 at the mean-field level, respectively. 
Using the above improved values of fi/Ep and A/Ep in the expression ([BT)]) . one obtains q c /kp ~ 0.349 
which is improved with respect to the mean-field value but still larger than the experimental one. On the 
other hand, the theoretical value 0.67 of the exponent £ obtained from Fig l26l at unitarity turns out to be 
quite close to the experimental value 0.66 (while the corresponding theoretical value for the prefactor T is 
2.5). Note that the comparison between theory and experiment proves better for the exponent C than for the 
prefactor T and the value of ^ andau in Eg. ([72"]) . This reminds us of the known fact about scaling relations, 
for which only the critical exponents and not the prefactors are considered to be universal. 
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It should be mentioned in this context that the true "critical" exponent Q c of a relation of the kind ([72")) 
where Vq — > + strictly (meaning that the lower bound 0.005 *~ Vo/Ep considered for obtaining the exponent 
Q reported in Figl26l is now removed) can be obtained analytically at unitarity as well as in the BCS and 
BEC limits by means of a local-density argument similar to that considered at the end of sub-section 4g. 
Specifically, one obtains ( c = 1 in the BCS limit as well as at unitarity, and £ c = 1/2 in the BEC limit. 
Comparison with Fig |26l then shows that £ = Cc in the BCS and BEC limit, while at unitarity the two 
exponents differ from each other. This difference is due to the fact that at unitarity the critical exponent 
C c can actually be revealed only for extremely small values Vo/Ep ~ 0.001. Approaching the BCS limit, 
on the other hand, the range of Vo/Ep where the critical exponent £ c holds becomes progressively wider, 
extending beyond the lower bound Vo/Ep — 0.005 considered for Fig(2H Along these lines we expect the 
value Cc = 1 to be the imprint of the fermionic character of the system, which is restricted to a progressively 
smaller range of Vo/Ep upon approaching the BEC limit, when the bosonic value ( c = 1/2 takes effectively 
over for all practical purposes. 

4g. Landau-type criterion for Andreev bound states 

In sub-section 4e we have emphasized the role played by the Andreev bound states on the BCS side 
of the crossover up to unitarity, in enabling the supercurrent to flow between the two sides of a barrier. 
These bound states, localized about the position of the barrier, act as a "bridge" through which most of the 
local transport occurs. The energy levels corresponding to these bound states lie below the threshold of the 
continuum levels. A typical Andreev level is shown in Fig |27[ with reference to the continuum branch of 
FiglUa) identified by the left-hand side of Eq.®. 




Figure 27: An Andreev bound state with ky = (red line) lies below the threshold of the continuum (blue 
line) corresponding to FigUJa). Here, a supercurrent with q/kp = 0.05 flows through a rectangular barrier 
of width Lkp = 4.0 and height Vo/Ep = 0.4 when the coupling is (kpap)^ 1 = —1. 



In addition, in sub-section 4f we have discussed how the self-consistent solutions of the BdG equations 
(which we have obtained across the BCS-BEC crossover) are sensitive to the instability associated with 
the Landau criterion for superfluidity when the barrier is allowed to vanish (Vq — > + ). In particular, on 
the BCS side of the crossover up to unitarity, we have recovered the onset of the instability due to the 
spontaneous generation of pair-breaking excitations, whose energy becomes negative as soon as the wave 
vector q associated with the superfluid current exceeds the critical value q c given by Eq.(|SD)). This limiting 
situation corresponds to the vanishing of the threshold of the continuum levels (in Fig(27J this threshold is 
represented by the bottom of the left dip of the curve). 
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Figure 28: Energy of the Andreev bound state (in units of Ep) with kii = vs q/kp when (kpap) 1 = —0.8, 
for rectangular barriers with given width Lkp = 4.0 and different heights. 
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Table 2: Ratio of the critical values of q, obtained from the maximum of the Josephson characteristics and 
the vanishing of the energy of the Andreev bound state with k|| = 0, for given barrier and coupling on the 
BCS side. 



In the presence of a finite barrier (with Vq ^ 0) when an Andreev bound state lies below the threshold of 
the continuum, upon increasing the value of q it is clearly the energy of this bound state to reach the zero 
value before the bottom of the continuum would do so. In this situation, it is natural to wonder whether an 
analog of the Landau criterion could be formulated even in the presence of a finite barrier. 

With reference to the BCS side of the crossover where the Andreev bound states are expected to play 
a crucial role, we thus attempt to formulate the "conjecture" that, in the presence of a finite barrier, the 
superfluid continues to flow until the energy of the lowest Andreev bound state (with k|| = 0) becomes 
negative. 

To verify this conjecture, we follow the evolution of the energy i?Andrccv of this bound state for a given 
barrier while increasing the value of the wave vector q of the supercurrent, and determine accordingly whether 
there exists a critical value g Androov at which i?Androcv vanishes (within some numerical error). To this end, 
in Fig |28l we report i?Andrcov vs q for several barriers (with the same width but different heights) for the 
coupling value (kpap)^ 1 = —0.8 on the weak-coupling side of the crossover, while in Fig[25] we do the 
same for barriers with the same height but different widths (and for the slightly different coupling value 
(kpap)^ 1 = —1.0). The corresponding values of g Andreev can be readily extracted from these plots. 

The values of g Andreov extracted from these plots are then compared in Table 2 with the critical values 
^josephson corresponding to the maximum of the associated Josephson characteristics for the same barrier. 
From this comparison we conclude that these two values of q c are indeed in quite good agreement with each 
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Figure 29: Energy of the Andreev bound state (in units of Ep) with kn = vs q/kp when (kpap) 1 = —1.0, 
for rectangular barriers with the same height Vo/Ep = 0.2 and different widths. 



other within some numerical error. This result then verifies explicitly our conjecture about the role played 
by the Andreev bound states for the collapse of the superfluid flow, at least for the couplings here considered 
on the BCS side of the crossover. 

This conjecture, however, starts failing as soon as the fermionic coupling becomes stronger and the 
unitary limit is approached. This can be seen from Figl30l where the energy i^Andreev is reported vs q/kp 
for a given barrier and several coupling values. It appears evident from the figure that i^Andreev is far from 
vanishing when the maximum sustainable current is approached, as soon as the coupling (kpap) -1 enters 
the crossover regime. This occurs in spite of the fact that the fraction of the current carried by the Andreev 
states becomes maximum at unitarity. 

Nevertheless, the point we like to emphasize here is that, for a given barrier, the coupling at which the 
energy of the Andreev bound state no longer approaches zero for the maximum sustainable current has 
already effectively verged toward the BEC side of the crossover. This can be clearly seen in Fig[3l] where 
the critical values q c corresponding to the four curves of Fig[3D]are plotted vs (kpap) -1 , and the curve that 
interpolates over these data points is compared with the two Landau limiting branches of Figl2"4l 



This comparison shows that the maximum of the curve corresponding to the finite barrier has moved 
toward weaker couplings with respect to the limiting Landau curve associated with a vanishing barrier. In this 
way, features of the BEC side of the crossover are now effectively felt already for couplings (kpap)" 1 < 0, so 
that the Andreev bound states (for which the internal structure of the Cooper pairs is an essential ingredient) 
loose the special role they have on the BCS side of the crossover. 

It is worth mentioning in this context that, for coupling values such that the critical current q c in the 
presence of a barrier is no longer determined by the occurrence of the Andreev bound states localized at the 
barrier, an alternative mechanisms is activated inside a wide enough barrier, whereby it is the local sound 
velocity to limit the amount of supercurrent that can flow across the barrier. Suppose, in fact, that the 
barrier is wide enough for the relation j = n(x)v(x) to occur between the local density n(x) and velocity 
v(x), in such a way that j does not depend on x. As a consequence, the local velocity v(x = 0) = j /n{x = 0) 
at the center x = of the barrier can be quite considerably increased over its bulk value j/no away from 
the barrier, owing to the local depression (n(x = 0) < uq) of the density at the center of the barrier, in 
such a way that v(x — 0) may as well exceed the local value s(n(x = 0)) of the sound velocity. Under 
these circumstances, the Landau criterion is activated locally and the superfluid current ceases to flow. The 
following simple relation 

q c _ n(x = 0) s(n(x = 0)) 
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Figure 30: Energy of the Andreev bound state (in units of Ep) with k|| = vs q/kp, for a given rectangular 
barrier with Lkp = 5.30 and Vo/Ep — 0.20, and several coupling values across unitarity. The percentage of 
the current carried at the upper value of q by the whole branch of Andreev states at the barrier center is 
also reported. 



is then found to hold between the critical wave vector for a finite (q c ) and a vanishing (gLandau-j bamgr (in 
the BEC limit, the right-hand side of Eg. ([73") ) reduces to = 0)| 3 /|$ | 3 where $(a;) is the solution of the 
GP equation). 

To show the (approximate) validity of this local criterion, we report in Fig[32] the ratio q c /q^ andliu for 
different coupling values (kpap) -1 on the BEC side of the crossover for a wide enough barrier, as obtained 
alternatively from the maximum of the Josephson characteristics and from the right-hand side of Eq. (|73|) , 
where the local density n(x — 0) results from the solution of the BdG equations © and the local sound 
velocity s(n(x = 0)) is obtained from Ref.[JS] for the given coupling and value of the local density. Note 
that the two values of q c are comparable to each other, even though the agreement becomes worse when 
approaching the strong-coupling limit. This is consistent with the observation we made in sub-section 4b, 
concerning the dominance of quantum tunneling when approaching the BEC side for given barrier, since 
the barrier height Vo eventually becomes larger than the bosonic chemical potential [is- In addition, when 
approaching the BEC side the healing length (cf. the quantity £ p hasc of Fig[7]) becomes larger than the 
barrier width, thus further invalidating any local approximation. 

We mention in this context that an analogous local Landau criterion was recently considered in Ref.[5Tj 
within an hydrodynamic approach based on a local-density approximation (see also Ref.[82 ). 



5 Stability of the self-consistent solutions 

In the previous Section, we have reported a number of numerical results that were obtained by solving 
the zero-temperature BdG equations throughout the BCS-BEC crossover until self-consistency was attained. 
Even though our numerical procedure has proven capable of determining these self-consistent solutions rather 
unambiguously, the question naturally arises about the stability of these solutions when fluctuations beyond 
mean field would be introduced. 

To a certain extent, this question can be addressed still remaining within our mean-field treatment, by 
exploiting the knowledge of the non-self-consistent solutions of the BdG equations about the self-consistent 
one. In a related fashion, one may also regard the two approaches (namely, of finding the mean-field solution 
via a self-consistent calculation as we did so far or via an energy minimization as we are going to do) as 
complementing each other along the lines of standard alternative mean- field approximations 83J . 

In the following, we shall first provide some general considerations concerning the expression for the 
energy related to the (not necessarily self-consistent) solutions of the BdG equations in the presence of a 
current, before embarking in its explicit numerical calculation for the cases of interest. The BEC limit, for 
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Figure 31: The curve that interpolates the critical values q c vs (kpap) -1 corresponding to the four cases of 
FigEOl is compared with the two Landau limiting branches of FiglMl 

which an independent knowledge of the results obtained by the GP equation can be exploited, will again be 
especially emphasized. In particular, we shall verify that the critical value q c obtained (like in Fig |25p for 
given coupling and barrier corresponds also to the maximum value of the superfluid current for which the 
self-consistent solutions of the BdG equations are energetically stable. 

5a. General considerations about the energy stability 

The quantity that needs to be estimated for the purpose is the total energy £ of the system, as it changes 
when non-self-consistent solutions about the self-consistent one are considered. Specifically, we have to 
compare the total energy of the system in the presence or in the absence of a barrier when the respective 
solutions are subject to the same "macroscopic boundary conditions" of given superfluid current j and overall 
phase difference A<p. The reason for taking this energy difference is that, even to a strictly homogeneous 
solution in the presence of the superfluid current j, there corresponds an energy cost needed to set up the 
current and the finite phase difference A(f> = 4>{xb) — 4>{xa) between two distant spatial points xa and Xb 
with (xb — xa) = D £ p hase |84j . [Note that, for notational convenience, in the present Section the phase 
difference A<p corresponds to the quantity 8<f>/2 of Section 4.] This energy cost should be appropriately 
subtracted off from the corresponding value of the total energy in the presence of a barrier (while keeping 
the same values j for the superfluid current and A(f> for the overall phase difference that accumulates accross 
the barrier), in order to probe the stability of the solution when a barrier is introduced into the system (say, 
in the middle of the interval D). 

By this procedure, we mean that the reference solution of the homogeneous case in the absence of a 
barrier (nb) should not only carry the same value j = qno/m of the superfluid current (no being the bulk 
density) but should also contain an infinitesimal local phase gradient (f>' n ^(x), in such a way that the total 
accumulated phase difference 4> n b(xB) — 4>nb(%A) has the same value A<f> as for the corresponding solution in 
the presence of the barrier. The only difference is that the latter solution develops the total phase difference 
A<j) sharply in a rather narrow spatial range of the order £ p hase about the position of the barrier, as we have 
explicitly verified by our numerical calculations. 

In addition, since the value of the chemical potential \i is kept uniform throughout the system, rather 
than the total energy £ we should actually consider the value of the grand-canonical energy £ — fj,N in the 
frame where the barrier is at rest, and look for the minimum of £ — fiN to determine the most probable 
configuration of the system 85J. 

To obtain the value of the total energy £, one could be tempted to calculate it by using the following 
standard expression of many-particle theory [86) 

£ = fdv hm i J2 e ^"" + V ^ + lLJ ") 0n( r > r >n) (74) 
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Figure 32: Ratio q c /q^ ndan vs (k F a F ) _1 on the BEC side of the crossover for a rectangular barrier of width 
Lkp = 2.65 and height Vq/E f = 0.10, obtained from the solution of the BdG equations (full line) and the 
expression (|T3")) (dashed line). 



which holds in the present form for a temperature T below the critical temperature. In this expression, 
/3 = l/(fcsT), u) n = (2n + 1)tt//3 (n integer) a fermionic Matsubara frequency, r\ a positive infinitesimal, and 
Qn the "normal" component of the single-particle Green's function in the broken-symmetry phase which 
admits the following eigenfunction expansion in terms of the solutions of the BdG equations [50] : 

Sn(r,r>„) = £ MrT + Mr') \ 

n v ' 

The problem with the expression (j?4")) when using the solutions of the BdG equations in Ea. (jT5"|) . however, 
is that Ea. (|74p holds when the self- consistency condition for these solutions has been achieved, while for the 
present purpose we need to know the value of £ also away from self-consistency. This important point results 
particularly evident when one attempts to derive the strong-coupling (BEC) limit of Eq. (|74]) . which is seen 
to coincide with the expression resulting from the Gross-Pitaevskii theory [87] only when the self-consistency 
condition represented by the GP equation (jSTj) itself is applied. The correct derivation of the BEC limit for 
£ also away from self-consistency will be considered in detail below. 

To obtain the expression of the total energy that holds also away from self-consistency, we have to resort 
to the general procedure discussed in Ref. pQ and, besides the original (grand-canonical) Hamiltonian 

H -fiN = J2j ( r )H(r)* CT (r) - g J dr*}(r)*}(r)*i(r)* T (r) (76) 

expressed in terms of the field operator x f , cr (r) (a =t, J, being the spin component) with a contact interaction, 
we have to introduce the effective (grand-canonical) Hamiltonian: 

H cB -^N = /**t(r)K(r)*,(r) 

+ Jdr [A(r)*j(r)*{(r) + A(r)** x (r)* T (r) . (77) 

Here, A(r) is the same arbitrary function entering the BdG equations ([5]) before the self-consistency condition 
© is eventually achieved. We then evaluate the thermal average ((H — fJ^N))jj eS _ m at with respect to the 
effective Hamiltonian ([77]) . and obtain the following expression with the use of the BdG equations ([5]): 



{2X> o 



((H -fiN)) HeS ^ N = rfr 2Ve„ [\u v (r)\ 2 f F (e u )-\v v (r)\ 2 (l-f F (e v ))] 
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2A(r) ^ U „(r)*^(r)(l-2/ F ( e ,)) 
g^u v (T)*v v (T)(l-2f F ( i e v )) 

(78) 



X>*,(r)Mr)*(l-2/,M)j 



This expression is formally a functional of the function A(r) that appears in the BdG equations ([5]), a 
function which can a priori assume arbitrary profiles when the self-consistent condition ([5]) is not satisfied. 

In particular, in the zero-temperature limit (whereby /f(e) = for e > 0) which we restrict to in the 
present paper, Eq. (|75|) simplifies as follows: 



fdv J -2^ e„Mr)| 2 + 2 A(r) £ u v {vfv v {v) 
ff ^u,,(r)*t )v (r)5]v(r)v(rr| . (79) 



The above expression need be regularized owing to the presence of the contact potential with strength g. 
This can be done by assuming that the non-self-consistent solutions of the BdG equations share the same 
ultraviolet behavior that was identified in order to regularize the gap equation (|13[) . In the zero-temperature 
limit we then obtain in the place of the expression (|79p : 

(T=0) 



|dr|-2^ £ ,K(r)| 2 + |A(r)| 2 J j£ 



dk m 
f k 2 

|A (1 ,|^— - ) 



where use has been made of the relation ^ to eliminate g in favor of o,f- 

As we have already discussed at the beginning of this Section, what we have actually to compare are the 
values of the quantity £ — /j,N = ((H — p,N)) calculated in the presence and in the absence of a barrier, 
while keeping in the two cases the same values of the superfluid current j flowing in and out of the system 
and of the overall phase difference Acf>. We thus write with reference to Eq.i 



{£ - fiN) {£ - fiN) 




A 



no barrier 



(|A(x)| 2 - A 2 ) 



dk. m 
(2tt) 3 k2 



(lA^I 2 - A 2 )^! -jAJ, (81) 

where A is the area transverse to the current flow (that we here consider finite for convenience), is defined 
after Eq.([2]), and Ao is the bulk value of the gap parameter far away from the barrier. 

Note the presence in Eq. (|8ip of the term j A0, which represents quite generally the energy cost required 
to set up the superfluid current j and the phase difference A(f) |84j . In the present context, when superfluid 
fermions are described by the homogeneous BdG equations with no barrier, the origin of the term j A(f> can 
be understood as follows. 

For superfluid fermions flowing with no barrier, we can use the results (115l) - (|17|) with (j, = fio + q 2 /(2m), 
which imply that A q , it q (k), and u q (k) do not depend on q as long as the superflow is stable. As a 
consequence, the first term on the right-hand side of Eq. p^|) does not contribute to the right-hand side of 
Eq. ([8T)]) when specified to the homogeneous case with the eigenvalues e u replaced by E^. The only dependence 
of £ = (H) on q thus originates from the explicit dependence on q of the chemical potential on the left-hand 
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side of Eg. (|80[) . Writing the total particle number as N = nV, where n is the uniform fermionic density and 
V = AD is the volume of interest containing the system, we then obtain: 

o\q\ o|q| m 

Suppose now to vary |q| — > |q| + A by an infinitesimal amount A, in order to generate the overall phase 
difference Acf> while mantaining unchanged the values j of the local current and /i of the chemical potential. 
Owing to Eg. ([52"]) . this variation induces the following increase in the total energy: 



where j = qn/m and 



■IB 



Ac/) = XD = I dx\= f B dx d ^ {x) . (84) 
Here, 4> n b(x) represents the additional phase entering the fermionic wave function [cf. Eq. fTO)) ] 

v„(r) -> u q (k) exp {i(k - q) • r - i(p n b{x)} (85) 



which is needed to produce the overall phase difference Atfi even with no barrier. The result (|83[) justifies 
eventually the presence of the last term on the right-hand side of Eq. (|8Tj) . 

The expressions (1751) and ([50)1 are written in terms of fermionic quantities and are accordingly able to span 
the whole BCS-BEC crossover from the weak- to the strong-coupling limits. However, the same conceptual 
and practical instances that have led us in sub-section 2d to consider the strong-coupling (BEC) limit of 
the BdG equations, by deriving from them the GP equation ([STjl for composite bosons, manifest here for 
the calculation of the total energy which we would like to express directly in terms of the condensate wave 
function $(r) of composite bosons [cf. Eq.([5"2])], We then need to extend the results of Ref.[5"0] in order to 
provide a formal derivation of the expression of the total energy in terms of composite bosons, starting from 
the original expression in terms of the constituent fermions. We include here for completeness this missing 
derivation. 

To this end, it is convenient to express from the outset the grand-canonical energy (I78|) in terms of 
the "normal" (Qu) and "anomalous" (G21) single-particle Green's functions that are compatible with the 
solutions of the BdG equations for an arbitrary choice of the function A(r) (this is required by our demand 
of studying the stability of the self-consistent solutions of the BdG equations). We have already reported in 
Eg. ([75)) the form of the "normal" component On, which holds even when the eigenfunctions {u v (r), v„(r)} 
and the associated eigenvalues {e„} are solutions of the BdG equations ([5]) with an arbitrary A(r). The 
corresponding form of the "anomalous" component Q21 reads: 

Mr,r>n) = E h ir)MrT - Mr) *^ (r/) ) ■ (86) 

The expressions ([75|) and ([56]) allow us to rewrite the grand-canonical average ([75]) of interest in the alternative 
form: 

((H - nN^H^-^N = 2 fdr lim i V e^'«(r)fi u (r ) r / ;w B ) (87) 

To obtain the limiting form of this expression in the BEC regime, we resort to the two integral equations 
that couple Qu and Q21 [50j . namely, 

e?ii(r,r';o; n ) = G {r,r';u n ) + fdr" Q (r, r'>„) A(r") G 21 (r", r'; co n ) 

g 21 {v,v'-io n ) = - fdr"g (r" ) r;-uj n )A(r")*g u (T",r';uj n ) (88) 
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where the non-interacting Green's function Qq satisfies the following equation 

H(t)] g (r,r';uj n ) = 5(r - r') . 



(89) 



We then iterate these integral equations, to expand Qn and C/21 to the forth and third order in A, respectively. 
After a long but straightforward calculation and taking into account the results of Rcf . 50J , we obtain for 
the first term on the right-hand side of Eq. (|87[) the approximate expression: 



dr 
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8tt 
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8tt w 4m w 



(Mb - 2V(r)) 
6ttclf 



|A(r)| 2 



2 \ 2 
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A(r) 



(90) 



where /is is the chemical potential for composite bosons defined after Eq. (|5ip . In a similar fashion we obtain: 



'£ 2 i(r,r; w») 



m 2 ap 



8tt 



m 2 aF V 2 
87T 4m 
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( MB - 2F(r)) 
47rai? 



m 2 a,F 
8tt 



A(r)* 
2 



|A(r)| 2 A(r)* 



which yields for the second term on the right-hand side of Eq. (187|) the approximate expression: 




' ? ' £ 2 i(r,r;uv) 



m 2 aF . / n V 

^T A(r) 4^ A(r) 



A(r)| 4 



(91) 



(92) 



where only terms which survive the regularization when g — > have been retained. Upon combining the two 
results (|9TJ|) and (f9"2"| . recalling the relation (|52"j) between A(r) and the condensate wave function <&(r), and 
making explicit the relation 

'drl^r)! 2 



(93) 



that holds in the BEC limit, we obtain eventually the desired result for the total energy in this limit: 
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£ = (H) 



H,, 



-fiN 



n r ( v 2 

- e °W dr {-* (r) 2^* 

^(r)|$(r)| 2 + i^|$(r)| 4 

2 m_e 



(r) 



(94) 



Here, m^ = 2m is the mass of the composite bosons, Vb(i") = 2 V(r) is the external potential felt by the 
composite bosons, and = 2a^ is the scattering length associated with their residual mutual interaction at 
the present level of approximation [50) . Apart from the term —6qN/2 which accounts for the energy required 
to form N/2 fermion pairs each with binding energy eoj the above expression coincides with the energy of 
a non- uniform bosonic condensate within the Gross-Pitaevskii theory (cf., e.g., Ref.[57]). This result proves 
that the fermionic BdG equations are able to recover a sensible description of the system of composite bosons 
that form in the strong-coupling (BEC) limit, also as far as their total energy is concerned. 

As already remarked, the results obtained in terms of the GP theory not only constitute an important 
benchmark for analogous results obtained from the BdG equations in the strong-coupling (BEC) limit, but 
may also independently serve to suggest methodologies and approximations to be adopted also for the BdG 
equations. To this end, it is convenient to write the analogue of the expression (|81[) directly in terms of 
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(composite) bosons. With reference to the form (|118[) of the bosonic wave function (cf. Appendix C), we 
obtain: 
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(95) 



where Nb = N/2, gs — '{iras/niB, and /ib = Qb^b + a B/(^ m B) at the present level of approximation. 
Evaluation of the bosonic expression (|95[) should thus be fully equivalent to evaluation of the fermionic 
expression ([5T]) once the strong-coupling (BEC) limit of the latter is considered. 

The form ([9"5"]) coincides with the expressions (18) and (19) reported in Ref . [55] . where the same boundary 
condition appropriate to a constant flow was adopted. This boundary condition is enforced by the last term 
on the right-hand side of Eg. ([93)1 . which accounts for the term jA(f> [53] in the present context and whose 
origin can be readily understood as follows. For the bosonic condensate, we write in the absence of the 
barrier [cf. Eq. (HUD! : 

j(qBX + 4>nb( X )) 



$nb(x) = e l{qBX+,p ^ x » (96) 

where ub is the uniform bosonic density and the phase </>^ b (a;) = Agi (with \b <SC <7b) is needed to account 
for the overall phase difference A0. In this case, the bosonic part of the total energy ([M]) reduces to: 
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(97) 



The last term on the right-hand side does not depend on the size D of the system since it can be written as 
the product j A(j) [84] . where j = qB tib /mB and 



dx 



dx 



dx Xf 



XbD 



(98) 



5b. Numerical analysis of the energy stability 

A numerical analysis of the bosonic expression (j5S"]) to test the stability of the solutions of the GP equation 
for point-like bosons was already performed in Ref. [88] . albeit for a delta- like barrier only. Here, we begin 
our numerical analysis by extending the results of Ref. 88J to more general barriers of finite width, and then 
proceed by evaluating the corresponding fermionic expression ([5T]) so as to consider the whole BCS-BEC 
crossover. 

As an example, we then consider the solution of the GP equation ([ST]) in the presence of a Gaussian 
barrier of width cr/£ p hasc = 3.0 and height Vo/Ep = 0.7 (in fermionic terms, this corresponds to the coupling 
value (/iFflf) -1 = +3.0 such that /cF^phasc = 3 v / 7r/4). The resulting Josephson characteristic is shown in 
Fig l33f a). whose maximum occurs for 54> c /ir = 0.3 and <? c £phase = 0.345. Note that in the figure we have 
represented by a full line the "left" branch of the Josephson characteristic when < 5<fi < S(f> c and by a 
broken line the "right" branch of the Josephson characteristic when Scj) c < S(f> < tt, thus anticipating a 
relevant difference between the two parts of the Josephson characteristic. 



Such a difference results evident in Fig l33T b). where the quantity A£ identified by the right-hand side of 
Eg. ([55]) is reported for several values of q that are consistent with Fig|551a), by perturbing the self-consistent 
solution of the GP equation ([5T]) while keeping the value of the superfluid current uniform in space at the 
same value j of the self-consistent solution. Specifically, the perturbed form of the condensate wave function 
about the self-consistent (sc) solution of the GP equation (|5ip has been chosen as follows. Its amplitude is 
varied locally according to the expression 

1*0,01 = l<M,)| + r |<M* = 0)1 ' (99) 
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Figure 33: (a) Josephson characteristic q vs 8<j> obtained from the GP equation with a Gaussian barrier 
of width cr/£phasc — 3.0 and height Vq/Ef — 0.7 (the two branches, respectively to the left and right 
of the maximum, are evidenced by a different drawing); (b) Ratio of the energy differences given by the 
expression (|95|) . calculated respectively for the non-self-consistent (AS) and self-consistent (AS SC ) solutions 
of the GP equation, vs the variable r defined in Eq. ([M|) . The different curves correspond to the values 
g£ phase = (0.1, 0.2, 0.3, 0.345, 0.1, 0.2, 0.3) from top to bottom. 



while its phase 4>b(x) adjusts itself locally in order to keep the current uniform. Here, $o = y n o/2 is 

the bulk value of |$ sc (cc)| and r is a dimensionless parameter with a typical range < |r| ~ 1. In this 
way, the perturbed form is specified only in terms of the value of the varied solution at the center of the 
barrier (x = 0). This procedure, in turn, defines the independent variable r = (A|$|/|$ sc |) x= o = (\&(x = 
0)| — | &sc(x = 0)|)/|$sc(£ = 0)| that controls the energy variations reported in Figl337b). 

Other choices of variations about the self-consistent solution might as well have been considered. For 
instance, we could have kept the value of the self-consistent solution at x = unmodified and varied the 
spatial extension of the self-consistent solution, still preserving the value of the supercurrent uniform in space. 
[We found that the requirement of the current being uniform is essential for a correct identification of the 
energy stability about the self-consistent solution.] We again emphasize that the requirement of the current 
being uniform for all possible variations about the self-consistent solution necessarily forces the value 5<p of 
the total phase accumulated by the non-self-consistent solution across the barrier to change with respect 
to the value of the self-consistent solution. Correspondingly, this change need be considered also for the 
quantity A<fi — 5(f>/2 appearing on the right-hand side of Ea. (f§5")) and associated with the solution in the 
absence of the barrier. 

The curves AS(t) reported in Fig |3"37 b) correspond to several values of q associated with the self- 
consistent solution of Fig |3"37 a). We note, in particular, that the full curves in Fig |33f b) show a minimum for 
t = (which is associated with the self-consistent solution) for values of q such that Scf) < 8<fi c in Fig |3"37 a). 
On the other hand, the dashed curves show a maximum for r = for values of q such that S(j) c < 6<fr in 
Figj33fa). Correspondingly, we conclude that the left part of the Josephson characteristics is stable against 
fluctuations (at least, for the class of fluctuations here considered), while the right part of the Josephson 
characteristics is unstable. This conclusion is in agreement with the results of Ref . [55] for a delta-like barrier. 

In order to fine tune the actual value of q past which the self-consistent solution of the GP equation 
becomes unstable, we report in Fig l34l the behavior vs q of the coefficient ci of the polynomial expansion 

= 1 + c x r + c 2 r 2 + c 3 r 3 + c 4 r 4 + c 5 r 5 (100) 

which fits best the curves of Fig l3"37 b). [In general, the need to include terms up to the fifth order in r stems 
from the analysis of the analogous quantity in the BCS limit, but we leave it here also in the BEC limit for 
shorting up the discussion.] Results for two different Gaussian barriers are reported. 



One sees from these plots that c 2 remains positive as long as 6(j> < <5</> c , while it turns negative when 5<j) c < 5<f), 
thus confirming our expectation. [The feature that ci returns slightly being positive for 5<fi ~ n is interpreted 
as due to our numerical error, which is estimated from this deviation to be less than 3%.] 
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Figure 34: Coefficient ci of the polynomial expansion (|100|) vs <?£ p hase corresponding to the set of curves of 
Figl3"3lb). Results are reported for two different Gaussian barriers with Vo/Ep = 0.7 and: (a) cr/£ p hasc = 1-0; 

(b) cr/£ p hasc = 3.0. 



Once the stability of the self-consistent solutions of the bosonic GP equation (fSTj) has been established in 
the presence of a supercurrent for barriers of finite width, we pass to consider the self-consistent solutions of 
the fermionic BdG equations (O for arbitrary couplings across the BCS-BEC crossover and their variations 
from self-consistency. To this end, we consider variations of the local gap parameter A(x) still of the form 
with A(x) now replacing ^(x), which in turn leads to modified values of the wave functions u„(x) and 
V v {x) (with the phase 2<f){x) of the local gap parameter being also modified in such a way that the value of 
the supercurrent is kept uniform). 




Figure 35: The quantity A£/ A£ sc vs r obtained (dashed lines) from the solution of the bosonic GP equation 
as in Fig |33f bh is compared with the corresponding quantity obtained (full lines) from the solution of the 
fermionic BdG equations for the coupling (kpap)^ 1 = +3.0 and the values q/kp = (0.01,0.04,0.06) from 
top to bottom. 



Following a procedure that we have adopted several times in this paper, we begin by exploring the BCS- 
BEC crossover for the quantity A£{t) by considering again a coupling value ((/cfOf) 1 = +3.0) deep in the 
BEC region, in such a way that a direct comparison can be established with the results obtained previously 
via the GP equation. 

We thus report in Fig[3S]the results for A£(t) using a Gaussian barrier of width akp — 5.0 and height 
Vo/Ep = 0.03, obtained for different values of the supercurrent and using the GP equation (dashed lines) 
as well as the BdG equations (full lines). The good agreement resulting from these two sets of curves gives 
us confidence in the accuracy of the solution of the BdG equations even as far as the quantity A£(t) is 
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Figure 36: The quantity A£ / A£ sc is plotted vs r for a Gaussian barrier of width akp = 4.0 and height 
Vo/Ep = 0.1 and three values of the coupling (kpap) -1 across the BCS-BEC crossover: (a) +1.0; (b) 
0.0; (c) -1.0. For each coupling, the results corresponding to three values of q/kp are reported (from top to 
bottom): (a) 0.01, 0.04, 0.06 (with q c /k F = 0.061); (b) 0.01, 0.15, 0.20 (with q c /k F = 0.21); (c) 0.01, 0.05, 0.07 
(with q c /k F = 0.076). 



concerned. Here and in the following, only values of q such that 5<fi < 8<j) c corresponding to the stable 
solutions will be considered. 

Accordingly, we report in FigJJSl the curves AE(t) for the three characteristic couplings (kpO-F) 1 = 
(+1.0,0.0,-1.0) across the BCS-BEC crossover, using a Gaussian barrier of width akp = 4.0 and height 
Vq/E f = 0.1. Curves corresponding to three different values of q are reported for each coupling. In all 
cases, we have verified that the energy stability of the solutions is lost as soon as q approaches its critical 
value q c corresponding to the top of the Josephson characteristic for the given barrier. More precisely, to 
extract the accurate values of q at which the instability in A£(r) begins to manifest, we have also studied 
the g-evolution of the coefficient C2 in a polynomial fit of the type (jlOOp . 

The above analysis can be repeated for progressively smaller barrier heights, up to the point that the 
conditions occur for the Landau criterion discussed in sub-section 4f to emerge. To this end, we plot in 
FigJ37]the behavior of the coefficient C2 of the polynomial expression (|100[) vs q at unitarity for progressively 
decreasing values of Vo/Ep, while keeping the width akp of the Gaussian barrier fixed. We see from this 
figure that, when the critical value q c (at which the coefficient C2 changes sign) approaches the Landau value 
<?Landau as soon as Vo/Ep — > + for the given coupling, the curves A£(r) become indeed progressively more 
unstable (this process becomes actually quicker as soon as q^ andau j g a pp r0 ached) . 




Figure 37: Coefficient C2 of the polynomial expansion (|100p vs the wave vector q of the supercurrent at 
unitarity for a Gaussian barrier of width akp = 4.0 and decreasing values of the height Vo/Ep: 0.1 (red 
squares); 0.01 (green triangles); 0.001 (blue circles). [The data have been multiplied by the indicated factors 
in order to make a common drawing of the three curves.] The critical value q c /kp = 0.380 at which C2 
changes sign in the limit Vo/Ep — > + compares favorably with the corresponding value g^ andau /fci? = 0.397 
resulting from the Landau criterion at unitarity. 
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These results conclude our analysis about the energy stability of the self-consistent solutions of the BdG 
equations in the presence of a barrier across which a supercurrent flows. This analysis, in turn, reinforces 
the validity of these mean-field solutions, to the extent that they are thus expected to somewhat survive the 
inclusion of fluctuations beyond mean field. In addition, the same analysis, when extrapolated to vanishingly 
small barriers, has also been proved to be consistent with the insurgence of the Landau instability in our 
numerical solutions of the BdG equations. This is a priori a nontrivial result, which then reinforces the 
validity of both analyses we have made about the energy and Landau instabilities. 



6 Concluding remarks 

In this paper, we have presented the results of a systematic numerical solution of the Bogoliubov-de Gennes 
equations at zero temperature for superfluid fcrmions with an arbitrary mutual interaction, adopting a spatial 
geometry (the presence of a barrier of finite width) and a boundary condition (the presence of a steady 
superfluid flow) which are characteristic of the Josephson and related effects. Reaching full self-consistency 
in the calculations has proved essential to obtain meaningful results, especially when approaching the BEC 
side of the crossover. 

This work represents the first systematic study of the Josephson and related effects throughout the BCS- 
BEC crossover, a topic which is actively studied at present both experimentally and theoretically. Extending 
the study of the Josephson effect with continuity from the BCS (weak-coupling) limit, where it has been 
traditionally confined, across the unitary (intermediate-coupling) region and up to the BEC (strong-coupling) 
limit, where coherence effects based on (composite) bosons eventually emerge, has allowed us to cast the 
physics of the Josephson effect under a more general and sometimes different perspective. Specific questions 
can be focused on in this context, which would not even emerge otherwise. We refer, in particular, to 
our findings about the role played by the Landau criterion for the lost of superfluidity when the relevant 
elementary excitations are generated in the limit of a vanishing barrier, and about the role played by the 
Andreev bound states in the presence of a finite barrier. 

Our emphasis on the achievement of self-consistency in the solutions of the BdG equations has made the 
nonlinear effects typical of a collective flow to emerge through the self-consistent requirement itself, while 
at each step of self-consistency the BdG equations were solved by scattering methods which are standard in 
the solution of the linear Schrodingcr equation. 

Implementation of the computational scheme has required us to develop rather sophisticated numeri- 
cal procedures. Also for this reason, it was important to rely on independent "benchmarks" with which 
the outcomes of the numerical calculations could be confronted. In the BCS (weak-coupling) limit of the 
crossover, when the size of the Cooper pairs by far exceeds the barrier width, the requirement that all bar- 
riers should be assimilated to a Dirac-delta barrier has represented a stringent test on the calculation. In 
the BEC (strong-coupling) limit of the crossover, on the other hand, the results of the numerical solution of 
the fermionic Bogoliubov-de Gennes equations have been favorably confronted with those obtained by the 
independent solution of the Gross-Pitaevskii equation for the composite bosons that form in that limit. This 
joining of the physics of the fermionic BdG equations to that of the bosonic GP equation represents per se 
an important step forward toward a unified understanding of coherence effects related to the broken gauge 
symmetry, especially as far as their practical implementation is concerned. 

The results obtained in this paper could be extended along several directions, besides a straightforward 
application of the present approach to more complicated geometries. In this context, an extension of the 
present approach to finite temperature, while still spanning the fermionic coupling from the weak (BCS) to 
the strong (BEC) limits, would require one to consider (at least) pairing fluctuations beyond mean field. In 
this way one should be able to recover the correct physics about the Bose-Einstein transition temperature 
in the BEC limit, as well as the physics of the pseudo-gap above the critical temperature in the unitary 
region. To achieve this goal, a formulation of the theory in terms of spatially dependent single-particle 
Green's functions should be preferred, in an analogous fashion to what we did in sub-section 5a relatively to 
the expression for the total energy in the BEC limit. 

Consideration of the time dependence for the Josephson effect would also considerably broaden the impact 
and application range of the theoretical results. This is partly due to the experimental setup utilized with 
ultra-cold Fermi atoms, whereby the presence of localized traps makes it preferable the use of initial and/or 
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boundary conditions leading to the onset of time-dependent effects. In addition, a general study of non- 
equilibrium time-dependent phenomena, which are intrinsically related to the building up of coherence and 
the onset of decoherence in macroscopic quantum systems, would constitute a significant task that could 
focus on the importance that macroscopic quantum phenomena are expected to acquire in the planning of 
future technology. In this context, one may anticipate that the time-dependent versions of the BdG equations 
for fermions and the GP equation for bosons will serve as starting points for the study of time-dependent 
macroscopic quantum phenomena throughout the BCS-BEC crossover, along similar lines to those followed 
in the present paper for the study of time- independent phenomena. In this respect, the stationary results 
obtained in the present paper could be used as initial and/or final conditions for time-dependent effects. 
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8 Appendix A. A conditionally convergent integral 

In sub-section 2b we have discussed in detail the procedure for constructing a complete set of eigensolutions 
of the BdG equations ([5]) in the presence of a supercurrent, to be used in the gap equation (or self-consistency 
condition) ^ as well as in the number ((5]) and current © equations. Apart form the bound states which 
may exist below the continuum threshold, in each of the remaining five energy ranges specified in FigslDJa) 
and 01b) we have used the appropriate wave vectors k^ in the outermost spatial intervals to label the 
continuum eigensolutions subject to outgoing boundary conditions. 

Under these circumstances, we have found it appropriate to normalize the continuum eigensolutions in 
terms of the wave- vector components k = k x = k± and {k y , k z ) = kii, where _L and || identify the directions 
"parallel" and "orthogonal" to the surface that fixes the geometry of the three-dimensional barrier (recall that 
k± coincides, by definition, with kg of Ea. (|22p in the outermost intervals where Vi = 0). In the numerical 
solutions of the BdG equations the wave vector k|| appears as a parameter through the definition of the 
reduced chemical potential (|2ip . The scattering problem is thus implemented by fixing the value of k|| first 
and then solving for the values of k± which extend in principle from — oo to +oo. 

With these premises, it would appear natural to perform the integral over the three-dimensional wave 
vector k = (fej_,k||), by integrating first over k± and then over k||. Nonetheless, the opposite order of 
the integrations could, at least in principle, be chosen as well. We shall show in the following, however, 
that the order of the integrations matters considerably, with the result that the intergral turns out to be 
convergent only when the integration over k± is performed first. On physical grounds, this difference stems 
by the different role played by the presence of the barrier when one is considering the two alternative limits, 
namely, of keeping k|| fixed while letting |A:_i_| to increase without bound, or doing the opposite of keeping 
|fc_i_| fixed while letting k|| to increase without bound. In the first situation the presence of the barrier should 
become progressively irrelevant, while in the second one the barrier keeps influencing the scattering wave 
functions no matter how large k|| is. 

Let us refer specifically to the gap equation at zero temperature in its regularized version f| 14[) . Apart 
from the discrete sum over the bound states (if any) , the sum over the index v therein reduces to an integral 
over k = (fcj_,k|i), where in the integral over A:j_ (which extends from — oo to +oo) the integrand takes 
alternative forms depending on the energy range it is spanning. For instance, in range 4 of FiglJJa) when 
fl > 0, k± coincides with k^ in the left branch and with fcW in the right branch of that figure, and the 
wave functions take the form ([5o]) or (|55]l . in the order. A similar situation occurs in range 7 of FigH^b) 
when fi < 0. 

To evidence the effects of the barrier in the two alternative situations discussed above (that is, of either 
keeping k|| fixed while letting |fc^| increasing without bound or doing the opposite), we plot as an example 
in FigJ35] the (square magnitude of the) coefficient b of the wave function (|25[) in the central spatial interval 
(namely, at x = where the potential barrier is maximum), which corresponds to an electron-like excitation 
impinging from the left (that is, with a\ = 1 in the outermost left interval). We see from this figure that, 
when k|| is kept fixed and k± is allowed to grow without bound (dashed line), \b\ 2 becomes exponentially 
small, thus showing that the barrier has no influence on the wave functions in this limit. In the reverse 
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situation, when k± is kept fixed and |k||| is allowed to grow without bound (full line), \b\ 2 approaches a 
constant value that depends on k±_ and is in general not negligible. 

We consider first the convergence of the integral over k± in Eq. (|14p in the barrier region, for large values 
of \k±\ and given value of k||. In this limit, we have just concluded that the barrier does not influence 
the wave function appreciably. For instance, in the expression (|25p the coefficient a — > 1 while all other 
coefficients become negligible. It is then sufficient to consider the expressions ([251 ) for the coefficients ug 
and V£ entering the wave function TQ (x = 0; q, E), in which we rewrite: 



E + & 



(^+kf) 
2m 



+ V t - Mo + 



\ 



•(*? + kf) 



2m 



+ V t - Mo + A 2 



(101) 



[Note that far from the barrier where Ve = 0, Af = Aq, and kg = k±, this expression reduces to that for the 
homogeneous case with q = 0.] For large values of |fcj_| (and therefore of kg) we then obtain: 



Ai 
2e t 



4e? 



(102) 



where we have introduced the short-hand notation 
k\ + (fef — kj_) = k\ + 5|, we expand: 



(kj + k|)/(2m) + Vt - no. Setting further kj 
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(103) 



with k 2 = k\ + k 2 . Equation (|102[) thus becomes: 
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(104) 
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Figure 38: Magnitude square of the coefficient b of the wave function (f2l>)) in x — corresponding to setting 
ai = 1 in the outermost left interval. The two curves correspond to keeping fixed either k±_ — kp (full line) 
or |k||| = kp (dashed line) and taking the independent variable be |k||| or k± 7 respectively. For definiteness, 
we have considered the coupling value (kpap)^ 1 = —1.0 for which /i > 0, while fl can acquire both signs. 
The calculation has been done for a Gaussian barrier with akp = 3.5 and Vo/Ep = 0.2, and for a current 
with q/kp = 0.03. 
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Here, the first term on the right-hand side when inserted into Ea. lfH]) cancels the leading term of the 
second term on the right-hand side of that equation. In the numerator of the second term on the right-hand 
side of Eq. (|104[) . on the other hand, when |fcj_| is large enough we may neglect 8 2 /(2m) in comparison with 
Ve 0. In this limit, the leading term on the right-hand side of Eq. (Ti"4"l) reduces to: 

1 [°° „ , f°° „ A e V e 



dku ku / dk 



m 2 A t V t [°° 1 nn ,, 

_ (1 ° 5) 



ftjj — 2m /io 

where for convenience we have introduced in both integrals a common lower cutoff k c (such that k 2 > 2m/io 
when fiQ > 0) since we are here interested only in the ultraviolet convergence of the integrals. This result 
proves analytically that the double integral in Eq. (|14|l is definitely convergent when the integration over k± 
is performed first. We have also verified from the fully numerical calculation that the same k^ 2 power-law 
behavior of Ea. Q105]) results in the ultraviolet for the integrand of the last integral over ku in Eq.(|14p inside 
the barrier (outside the barrier, on the other hand, an even faster decay results since the subtraction of the 
homogeneous contribution in Eq. (|14[) is now more effective there). 

Let's, finally, consider the opposite case when k± is kept fixed and |kj|| is allowed to grow without 
bound. In this limit, we have argued that the effect of the barrier cannot altogether be neglected, since 
the scattering coefficients are in general non negligible for large values of k± where they approach constant 
values that depend on k± (cf. FigJ35]for \b\ 2 inside the barrier as an example). In this case, it is sufficient to 
demonstrate that the double integral over k± and ky in Eq. p4j) does not converge when the integration over 
ky is performed first, by considering that equation far from the barrier (for example, for x ~ — oo) as this 
was actually the most favorable situation for the convergence of the integral when the integration over k_L 
was performed first. We may also consider, for instance, the simple case of a rectangular barrier of width L 
and height Vb for which the scattering coefficients can be calculated analytically when kjj / (2m) is the largest 
energy scale in the problem. In the place of Eq. (|105p , we obtain the following expression for the leading 
term on the right-hand side of Eq. (T14"]) when x rs — oo: 

— - / dk^Reihik^ 2 ^} rife,, fey ° (106) 

W Jk c Jk c (2^ -Mo) 

(mV a ) 2 (l - cos2L-JE - 2mV ) 

Mk^)\ 2 = — ^ ? > , (107) 

2k j_ - 8k 2 ± (mV ) + (mVo) 2 (l - cos2L^/k 2 ± - 2mV ) 



where 



and we may assume k c > \/2mVo. It is clear from the expression p06p that the integral over fe|| therein 
diverges logarithmically in the ultraviolet. 

These results prove our initial statement that the order of the integrations over k± and ku in the gap 
equation (|14[) matters, owing to the different role played by the barrier when either k±_ or ku are allowed 
to grow first without bound. Consistently with these results, we have implemented the numerical solution 
of the BdG equations ([5]) by integrating first over kj_, as it would anyway appear more natural on physical 
grounds. 



9 Appendix B. An analytic integral for the total current with a 
Dirac-delta barrier 

We have seen in sub-section 2c while discussing the non- self- consistent solution of the BdG equations for a 
delta-like barrier, that the total current can be cast in the form [cf. Eq. (l4"7| ]: 

W) - [ *dfc||fcn = A °f n ^ [ kF dk l{ fc||-= Tk " (108) 
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where 

1 _ 1 _ 1 - k 2 
Tk a ~ i 7 72 - 1 1 m 2 z 2 - 1 1,2 • ( 109 ) 



Here, we have set k = ku/kp and 



T kF = = < 1 • (HO) 

1 + ^ 1 + z% - 



Simple manipulations then bring Eq. (|108p into the form: 

fefpAosin^ T kl 



47T , 2 



1 - T kF sin 2 8(t>/2 



1 n . i.2\ 



x / dkk ^ fc j (111) 

J0 r, 7p jTJ L T kp cos 2 30/2 , 2 



Identifying the two new quantities 

T fcF cos 2 (50/2 
1 - T kF sin 2 50/2 

and changing the integration variable from k to t = fc 2 , Ea. (|lll[) reduces eventually to the compact form: 



X = T kF and (US) 



fc| A sin 50 T kF f 1 (1 - t) 



J(5ct>) = ^-^ ^ = / di -L . (113) 

>\-T kr sin 2 SH2 J ° VT^XtVT^Yt 

In this expression one recognizes an integral representation of the Appell hypergeometric function of two 
variables [89] 

F 1 (a,P,P', T ,X,Y) = - -±>> / dt t a ~ l — -L_^_ _ (114) 



r(«)r( 7 - «) y (i-a)f (i-t r)/ 3 ' 

when a = 1, p = j3' = 1/2, and 7 = 3. In Eci. (|114[) . T(z) is the Euler's Gamma Function (such that 
r(n + 1) = n! for integer values of its argument) and by definition (a) n = a(a + l)(a + 2) • • • (a + n— 1) with 
ao = 1. Alternatively, the Appell hypergeometric function admits the double serie representation: 



^ ^ ( a ) n+m (P) m (p') n 



F x {a,p,p', r ,x,Y) = Y,22 "|"V" m ^ ; " A™F" (115) 
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which is absolutely convergent for \X\ < 1 and \Y\ < 1. In this way, the total current (| 1 1 3[) takes its final 
form: 

167r /l-T feF sin 2 50/2 V 2 2 ; 



with X and y defined in Eg . 1(112 



Recall that the expression (| 1 1 6|) holds for any value of Tfc F in the interval < T^ F < 1. In particular, in 
the limit of a weak barrier when T kp -> 1 Eqs. ((TT2"]) give X -> 1 and Y -»■ 1, so that Fi(l, 1/2, 1/2, 3; X -> 
1, y — > 1) = 2 and the expression (|116[) reduces to the result (|3U|) obtained in the text by a different method. 
Analogously, in the limit of a strong barrier when Tk F « .2^ 2 C 1 we may take X = and y = in the 
arguments of the hypergeometric function, so that Fi(l, 1/2, 1/2, 3; X = 0,Y = 0) = 1 and the expression 
(I116P reduces to the result (J49J) of the text. 

More generally, for intermediate values of Tt F the Josephson characteristics (|116[) can be obtained by 
computing numerically the hypergeometric function therein in terms of the series representation (|115[) . as 
shown in Fig |39l This plot shows the progessive evolution from J {5(f) oc sin 5(f) for a strong barrier (that 
is, small 7fc F ) to J {5(f)) oc sm5<f>/2 for a weak barrier (that is, large 7fc F ). Only the results obtained from 
Eci. (|116p for a strong barrier survive, however, the inclusion of self-consistency in the calculation, as discussed 
extensively in sub-section 2c of the text. 



10 Appendix C. Absence of the reflected wave in the Josephson 
tunneling 

In this Appendix, we analyze the flow of a Bose-Einstein condensate against a potential barrier at zero 
temperature as described by the Gross-Pitaevskii equation, by evidencing the analogies and differences with 
the scattering problem of a single particle against the same barrier as described by the ordinary Schrodinger 
equation. 

To this end, we begin by rewriting the Gross-Pitaevskii equation (|51|) entirely in terms of bosonic quan- 
tities as follows: 

V 2 

- $(r) + V B (r) *(r) + g B |$(r)| 2 *(r) = fi B *(r) (117) 

2m B 

where we have set g B = A-Ka B /m B . Were not for the nonlinear term proportional to g B in Eq. (|117p . the 
structure of this equation would be identical to that of the (time-independent) Schrodinger equation for a 
particle of mass m B subject to the potential V B {r) and with eigenvalue [i B . The nonlinear term, however, 
is responsible for the occurrence of phenomena which would not be possible according to the ordinary 
Schrodinger equation. 

In particular, in what follows we shall focus on the occurrence of a reflected wave when a current flow 
impinges against a potential barrier. While an ordinary scattering process against a fixed potential as 
described by the Schrodinger equation unavoidably results in a reflected wave (barring exceptional cases for 
special values of the energy), we shall argue that the flexibility introduced by the nonlinear term in the 
Gross-Pitaevskii equation may lead to a complete absence of the reflected wave. This is because the wave 
function $(r) is able to adapt itself locally in the vicinity of the potential barrier V B (r), with the result of 
shaping the induced potential proportional to g B in such a way to counteract the effect of the barrier as far 
as the the presence of the reflected wave is concerned. In this context, the absence of the reflected wave in 
the description of the condensate corresponds to the occurrence of the stationary Josephson effect of interest 
in this paper. It turns out that this effect can occur provided the current does not exceed a critical value, a 
condition which is also familiar for the Josephson effect. 

For a potential barrier with a slab geometry of the type considered in this paper, V B (r) depends only on 
the coordinate x orthogonal to the barrier and the solution of Eq. (|117D can be cast in the form: 

= V^bR ^ b{x) e lqBX (118) 

where n B {x) — |$(x)| 2 is the local bosonic density. [Note that, for a barrier with a slab geometry, the 
coordinates parallel to the barrier disappear from the solution of the GP equation when the current flow is 
taken orthogonal to the barrier. This marks a difference from the solutions of the fermionic BdG equations, 
for which the coordinates parallel to the barrier play a crucial role especially in the BEC limit where they 
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contribute in an essentially way to the formation of the composite bosons.] To the solution pi8p there 
corresponds the current 

m = ^( qB + . dig) 



TUB \ dx 

Far from the barrier, we expect n B {x) to recover its bulk value n° B and 4>b{x) to reach constant values 
(say, 4>b(x — ► — oo) = and 4> B (x — * +00) = 4>b), such that the current (|119[) is given by n B q B /mB 
(corresponding to a superfiuid of density n B flowing with velocity V = qs /tub)- To mantain this value 
everywhere, the phase 4>b{x) in Eq. (|119[) has to adjust itself to the local variation of the density ub(x). 
Under these conditions, the bosonic chemical potential in Ea. (|117p has the value \ib — 9sn B + 9s/(2ms). 
Note also here that in the absence of the barrier nothing would, in principle, prevent qs (and thus the 
current) to grow without bound owing to Galilean invariance. It is just the presence of the barrier, by 
breaking the translational symmetry of the system, to introduce an upper limit on the allowed values of qB- 

We thus consider the one-dimensional version of Eq. (|117[) 

' ' M> '' r) ' V B (x)$(x)+g B mx)\ 2 <S>{x) = iibHx) (120) 



2m b dx 2 
subject to the boundary condition [55] : 

= ^e lqBX e l,pB{x ^ ±oo) when x -» ±00 . (121) 

Suppose that we have found such a solution, which we label $ gB (2) for brevity. Recalling that (1b = 
9b i^B + 95/(2™^), we can manipulate Ea. (|120p and cast it in the form of a Schrodinger equation with 
eigenvalue q 2 B /{2ms)- 

1 d^ qB (x) + vM ^ {x) = ^_ ^ {x) (i22) 



2ms dx 2 B 2m b 

with the effective potential 



V ef£ (x) = V B (x) - g B n° B ( 1 - ) (123) 



B 

2 



where tib(x) — \Q qB (x)\ . Typically, |$ 9B (x)| (and consequently ub(x)) is lowered from its bulk value over 
a length scale £ p hase about the barrier, thereby producing an effective potential (1 1 23[) with the generic shape 
reported in Fig |40f a). The characteristic feature of this potential, which is relevant to our argument, is that 
it is repulsive in the central region and attractive outside. 

To make the solution of the Schrodinger equation (|122[) with a potential of this shape tractable, we 
further assimilate the potential to the simplified model form depicted in Fig l40f b). namely, 

f V if |z| < L/2 
V cS (x) = I -V x if L/2 <\x\<b (124) 
[ otherwise 

with Vq > and V\ > 0. For this potential, five spatial regions can be identified where the solutions of 
Eq. (|122p have the local form 

&t(x) = At e lkex + B e e~ lkex (I = 1, 2, • • • , 5) (125) 



with fei = k 5 = q B — \J2ttlbE, k 2 — k 4 — yj 2m B (E + V\ ) , and fc 3 = y/2m B (E — Vq). 



Requiring the continuity of the functions (|125p and of their derivatives at each interface, we can express 
the rightmost coefficients (A5, B$) in terms of the leftmost ones (Ai, B\) via the transfer matrix Y and write: 



A 5 \ f Y u Yj3 4i 
B 5 ) I Y 21 Y 22 J I Bi 



(126) 



with det Y = 1. For a wave impinging from the left we can set A\ = 1 and B$ = 0, yielding B\ = —Y 2 \jY 22 
and A5 = 1/Y 22 (the correct normalization A\ — > \Jn B will be imposed only at the end). In addition, 
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Figure 40: (a) Generic shape of the effective potential V c s(x) of Eq. (|123p when the external potential Vb(x) 
has the form of a rectangular barrier; (b) Simplified model potential (|124[) . showing the essential features of 
being repulsive in the central region and attractive outside. 



current conservation across the different regions requires |-Bi| 2 + l^sl 2 = 1, yielding 1 1 2 — |^2i| 2 = 1- In 
this specific example, the condition B\ — for the absence of the reflected wave (which we argued above 
to correspond to the occurrence of the Josephson tunneling) is then achieved when the matrix element Y21 
vanishes. 

Further progress can be obtained only by numerical calculations. To reduce the number of parameters in 
the calculation, and keep contact at the same time with the original problem, we take for Vb{x) a rectangular 
barrier of height Vq and width L and relate the values of Vq and V\ in Ea. (|124p to the original expression 
(I123|) by approximating therein Vq « Vq — gBn B {\ — n B (0)/n B ) and V\ « —g B n B {\ — n B (x)/n B ), where x 
is an average value of x in the attractive region of the effective potential. To simplify the treatment even 
further, we may also take ng(x) w ub(0) so that Vq « Vq — ag B n B and Vi ~ — ct g B n B with < a < 1. For 
fixed values of Vb, L, and gBn B , we can then solve for the scattering problem with a chosen set of values of a 
and <7b, thus determining the matrix elements of the transfer matrix (|126|) vs the parameter b of Fig l40f b). 
Quite generally, one finds that the equation = admits (infinitely many) solutions only when qs is 

smaller than an upper value q B , and that out of these solutions only the first two (that we label I and II) 
with < b^ < b^ 11 ' are acceptable, as they correspond to spatial profiles of the wavefunctions which are 
consistent with the assumed shape of the effective potential in Fig0O](a). 

As a specific example, we report the numerical results obtained for a barrier with Vq — 5gBn B and 
L = 1.5 £ p hasc where £ p hasc = 1/ V^™s5s^s- For such a strong barrier we expect us(0) ~ 0, so that we can 
take a = 1 in the parametrization of the effective potential, yielding Vq — 4,gBn B and V\ = —g B n° B . For 
a given value of gs, we first determine the two solutions b^\qB) and 6^ j '(?b) of the equation |£?i|(&) = 
discussed above, in correspondence to which we obtain two distinct values for the coefficient ^5 = 1/Y22- 
Writing further A§ — \As\ expji^s}, we eventually determine the two corresponding values <j^g and </)^ 7 ^ 
of the phase. 

By this procedure, we end up with the Josephson characteristic <7b(0_b) plotted in Fig|4"l7a). which can 
be nicely fitted by the expected relation qB{4>B) = Q B sin^s with q c B — 8.26 x 10 _3 C ph 1 asc - Note that for the 
same rectangular barrier, the numerical solution of the original Gross-Pitaevskii equation (|120p yields the 
same functional relation of q B vs <p B with the value q c B = 14.0 x lO -3 ^^,,. 



It is interesting to consider also the opposite case of a vanishingly small barrier, so that Vq = V\ = 
—agBn B and the effective potential of Fig |4llf b) reduces to a square well of depth |Vi| and width 2b. In 
this case, only three spatial regions need be considered for the scattering problem and one can readily solve 
analytically for the transfer matrix of Eq. (|126p . One obtains: 




or J q 2 B + 2m B \V 1 \ 
Y 21 = -\ QB - V ] sin(2&A / g | + 2m B |Vi|) (127) 

/q% + 2m B \V 1 \ qB 
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Figure 41: Josephson characteristic q B {4> B ) obtained from the solution of the Schrodinger-type equation 
(|122[) with an effective potential of the shape of Fig |4"07 b) and corresponding to a: (a) Strong barrier 
(whose parameters are specified in the text); (b) Vanishingly small barrier. In the latter case, the result of 
the approximate calculation (full line) is compared with the corresponding solution of the Gross-Pitaevskii 
equation (dashed line). 



which vanishes for b — — and b = b^ 11 ^ — n/ (2y q 2 B + 2ms|Vi|). Only the second solution is of interest 

here. In the present case, the value of the parameter a entering \ Vi \ cannot be determined a priori but need 
be determined self-consistently. To this end, it is sufficient to know the expression of the wave function in 
the central region, namely, 

2m B \Vi 



q 2 B + 2m B \V 1 \ 

whose average value in the interval —b^ 11 ^ < x < b^ 11 ^ equals: 



<S> 2 {x)\ 2 = 1 - 2 ■ ,'- . cos 2 {xJq 2 B + 2m B \V 1 \) , (128) 



q% + 2m B \Vi\ 



Recalling that by our definition a = 1 — (|$ 2 | 2 ), Ea. (|129[) yields eventually the result: 

1 Q% 1 



a = 



2 2m B g B n° B 



(130) 



where the assumption < a requires that q B < q% = \Zm B g B n g . Recalling that g B n B — m B s 2 , we thus 
obtain that the superfluid velocity V = q B /m B cannot exceed the sound velocity s of the original bosonic 
system, as required by the Landau criterion. 

Finally, from the matrix element Y\\ we can also determine the phase <p B of the non- vanishing coefficient 
in the rightmost region, and obtain for the Josephson characteristic the result: 



q 2 B + 2m B \V 1 \) V 9 S 

This expression is plotted in Fig HTT b) together with the corresponding result of the Gross-Pitaevskii equa- 
tion, which gives q B = q c B cos(</>b/2) for the same "second branch" of solutions. In spite of the crudeness 
of the model that we have adopted for the effective potential V c s{x), it is remarkable that our simplified 
solution captures the essential features of the exact solution, especially as far as the "second branch" is 
concerned (which always proves the most difficult one to be obtained). 
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11 Appendix D. Some useful manipulations for the numerical so- 
lution of the Gross-Pitaevskii equation 

It was stressed in sub-section 2d that one keystone of the present approach to the Josephson effect throughout 
the BCS-BEC crossover is the possibility of confronting the numerical calculations based on the fermionic 
BdG equations ([5]) in the strong-coupling (BEC) limit with the independent solution of the Gross-Pitaevskii 
equation (|51[) . provided some reasonable physical conditions are satisfied |50) . The conditions for the agree- 
ment between the two alternative (fermionic and bosonic) calculations were also discussed in sub-section 2d 
(cf., in particular, FigllOl therein), and were repeatedly explored in Section 4 while discussing the results for 
the Josephson and related effects with an SsS barrier. The relationship with the Gross-Pitaevskii equation 
was explored, in addition, in Appendix C for providing an alternative interpretation of the Josephson effect, 
whereby the role of self-consistency was emphasized as leading to a destructive interference for the wave 
reflected by the barrier. 

In this context, it is worth mentioning that the solution of the one-dimensional version (|120[) of the 
Gross-Pitaevskii equation which describes a steady flow past an obstacle was discussed with some length 
in the literature, even providing analytic results in some limits (cf. Ref.[88j and references quoted therein). 
As a matter of fact, these analytic results have been utilized in our calculations whenever possible. Yet, in 
most circumstances we had to solve numerically the equation (|120p in the presence of a barrier with a slab 
geometry and a finite current. To this end, some preliminary formal manipulations of this equation have 
proven useful to speed its numerical solution. We report them here for completeness. 

Recalling the boundary condition (112ip and the modified value of the chemical potential [i B = gBn- B + 
q B /(2m B ) in the presence of a current, we begin by dividing both members of Eq. (|120[ ) by g B ( n B ) 3 ^ 2 & n d 
introduce the notation £ p hasc = (4msgs n'g)" 1 / 2 , b(x) = Vb(x)/ gBn B , and q = <?b £ p hasc, as well as the 
rescaled variable y = x/£, phasc . Setti ng further $(y)/(n° B ) 1/2 = f(y) = f (y) exp{i(j) B (y)} exp{iqy} (where 
fo(y) is a real function - cf. Eo. (|118p ). the one-dimensional Gross-Pitaevskii equation (|120p then reduces to 
the form: 

2 /"(!/) + (l + 2q 2 - b(y)) f(y) - f Q (y) 2 f(y) = . (132) 

Here, the second derivative of f(y) can be expressed in terms of fo(y), </>s(y), and q, such that the real and 
imaginary parts of Eq. (|132p can be separated, in the order, as follows: 

2fS(y) - 2 f a (y) (q + cp' B (y)) 2 + f (y) (l + 2q 2 ~ b(y)) - f (yf = , , 

2/6(2/) (q + My)) + My) My) = o . 1 66} 

We note that the second of these equations corresponds to the condition of a uniform current. From the 
expression (|119D of the current, one writes, in fact, for this condition: 

fo(y) 2 (q + My)) = 9, (134) 

whose derivative yields the second of Eqs. (|133p . We can then conveniently replace the second of Eqs. (|133D 
by Ea.lUm 

Further manipulation can be done in order to remove the function <Pb(v) from the first of Eos. p33p . by 
making it to depend only on the function fo(y). Upon multiplying both sides of that equation by /o(y) 3 , 
the second term on its right-hand side becomes — 2/o(y) 4 (q + 4>' B (y)) 2 = — 2q 2 according to Eq. (|134p . In the 
place of Eqs. (|133p , we thus solve explicitly the two following equations: 



fo(y) 3 [2fZ(y) + My) (1 + 2? - b(y)) f (y) 3 ] = 2q 2 

k(y? (g + My)) = 9 ■ 

Once the first of these equations is solved for fo(y), the second one yields: 



(135) 



y , , (i - My 1 ?) 
My' 



•>B(y) =91 dy' f u (136) 



where we have assumed 4>B(y = — oo) = as usual. 

Finally, the order of the first of Eqs. (|135p can be lowered in the spatial regions where b(y) is constant. 
To this end, we multiply both sides of that equation by f$(y), such that it can be rewritten in the form: 

d f i( \2 , U + 2 q 2 - b (y)) d 2 1 d 4 _ 2 d 1 

^ /0(y) + 2 ^ ~ lTy M + q TyJW = ° • (13?) 
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In particular, outside a rectangular barrier b(y) — 0, so that Ea. (|137p reduces to the form: 

m 2 + (1+ 2 2e) My) 2 - \h{yf + q 2 j^ = \ + 2? , (138) 

where the value of the constant on the right-hand side has been obtained by considering either one of the 
limits y — > ±oo where fo(y — > ±oo) = 1 and fq{y — > ±oo) = 0. 

Once the numerical solution of the one-dimensional Gross-Pitaevskii equation in the presence of a current 
is obtained along these lines, comparison with the corresponding solution of the fermionic BdG equations 
in the strong-coupling (BEC) limit is eventually made by recalling the relations = 2m, nP B = no/2, 
qB = 2q, Vb(x) — 2V(x), and as — 2ap, as well as the rescaling (f52"j) between <f>(x) and A(x). 



12 Appendix E. Failure of the transfer-matrix method for solving 
the BdG equations 

When solving for an ordinary one-dimensional scattering problem described by the Schrodinger equation 
with a finite-range potential, it is often convenient to resort to the so-called transfer matrix whereby the 
amplitudes of the wave function at the farthest right of the potential are expressed in terms of those at 
its farthest left. In particular, when the profile of the potential is approximated by a step-like curve over a 
number of intervals as we did in sub-section 2b, we can solve the Schrodinger equation by elementary methods 
in each of these intervals and connect the solutions in contiguous intervals via the continuity conditions. In 
this way, the amplitudes of the wave function in the £-th interval can be expressed in terms of those of the 
(I — l)-th interval, and by recursion in terms of those of the interval with £ = 1 at the farthest left. Reaching 
eventually the interval with £ = M at the farthest right, one ends up with an operative definition of the 
transfer matrix. The advantage of this method is that one always deals with matrices of a small size (2x2 
for the ordinary Schrodinger equation); the disadvantage is that these matrices needs to be multiplied for a 
large (M — 1) number of times to obtain the required transfer matrix. 

When dealing with a one-dimensional scattering problem described by the BdG equations as in the 
present paper, it is then natural to ask whether it might be convenient to organize its solution in terms of 
the transfer matrix, which again considers the sequence of continuity conditions taken one at a time instead 
of grouping them altogether in a single equation as we did in sub-section 2b [cf. Eq. (|2"9")) ]. To obtain the 
transfer matrix, the continuity conditions J26|) at the point xg between the £-th and {£ + l)-th intervals have 
to be rewritten in the form 

Wt+i = M e+1 (x = xe)- 1 M e (x = x e ) W e (139) 
where We is the column vector (|2"7|> . such that eventually: 

W t=M = ( II M^'+iO* = a*) -1 M ?'( x = x i') ^ Wi=i (140) 

which defines the 4x4 transfer matrix for the BdG equations. The unknowns to be determined by solving 
for Ea. (|140p can again be identified by using "outgoing boundary conditions" as we did in sub-section 2b. 
For instance, when one selects an electron-like excitation impinging from the left with (<zi = 1, d\ = 0, 6 m = 
0, cm — 0), an electron-like excitation impinging from the right with (ai = 0, d\ = 0,6m = 1, Cm = 0), a 
hole-like excitation impinging from the left with (a\ —0,di = 1,6m = 0, cm = 0), or a hole-like excitation 
impinging from the right with [a% — 0, d\ = 0,6m = 0, cm = 1), in ah these cases the unknowns of the 
inhomogeneous system (|140j) are (6i, c\, ttjf, g?m)- When looking instead for bound-state solutions with (ax — 
0, g?i = 0, 6m = 0, cm = 0), the system (|140[) becomes homogeneous in the same unknowns c%, ojf, c2m)- 
All other coefficients in the intermediate regions with £ — (2, • • • , M — 1) can be determined, for instance, 
from (6i, ci) by applying the condition (|139p recursively from each interval to the next. Finally, the expedient 
introduced in Eqs. (f54"|) and ([55]) can also be used here for the same reasons. The expression (|139p is thus 
rewritten as follows: 

W e+ i = M £+1 {x = 0)- 1 Mi(x =x e - x^i) W e , (141) 

such that the values of the variable x entering the phase factors exp{ikj, x} in the expression (j28|) do not 
exceed the size of the largest interval. 
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By this procedure, one has traded the large [4(M — 1) x 4(M — 1)] size of the single matrix A in Eq.(J2U) 
with the [(M — 1) times] repeated multiplication in Ea. (|140p of a small [4 x 4] matrix. Although this trading 
may, in principle, seem favorable, in practice it turns out not to be the case because there occurs a severe 
limitation on the number n of iterations of the expression (|139[) that one can actually perform (and therefore 
on the maximun number M of zones one can rely on) . Following a similar line of reasoning to that introduced 

at the end of Section 3, we can take lm{k}/kF at most of the order 500, such that (500) 2 ™ ~ e 500 as required 
by standard numerical precision. This estimate yields n sa 40, which is indeed too small a number for a 
satisfactory determination of the gap profiles that has to be sustained in the cycles of self-consistency. The 
numerical attempts we have done along these lines have failed accordingly. For the above reasons, we have 
stuck with the method discussed at length in Section 2, which solves for all continuity conditions at once as 
evidenced by Eq. (|2"9)l . 



13 Appendix F. Failure of the tunneling Hamiltonian for describ- 
ing the Josephson effect in the BCS-BEC crossover 

The theory of electron tunneling in superconductors has historically been based on the use of the tunneling 
Hamiltonian [65] , following the original treatment by Josephson on coherent pair tunneling [90] . Only more 
recently the BdG equations have, in fact, been utilized to calculate physical quantities related to the Joseph- 
son effect [58 , especially as far as its connection with the Andreev-Saint-James bound states is concerned. 
In this Appendix, we wish to comment on the reasons why attempts to use the tunneling Hamiltonian are 
bound to fail when trying to describe the Josephson effect throughout the BCS-BEC crossover. 

Within the tunneling Hamiltonian formalism, the fermions at the left (L) and right (R) of a barrier are 
assumed to be independent from each other, in such a way that they are described by destruction Ck, CT and 
creation cj^ operators of wave vector k and spin component a (with k = k^ at the left and k = k# at 
the right of the barrier) which anticommute with each other. The barrier itself provides a coupling between 
these two otherwise independent systems, via the so-called tunneling Hamiltonian: 

H t = (*ktk« C L < r < W + *k i k*Ck J1 ,<r < w) ( 142 ) 

k L ,k H ,(T 

where tk L k R is the tunneling matrix element. A full account of the tunneling processes would then require 
one to specify the dependence of this matrix elements on its arguments in some detail. 

This problem is readily overcome in conventional superconductors, since the tunneling takes place over 
a narrow energy range of the order of the Debye energy about the Fermi surface. The fermions involved in 
tunneling have thus wave vectors near the Fermi wave vectors hp and kp on the two sides of the barrier, 
so that it is appropriate to consider the tunneling matrix elements in Ea. (|142p as being approximately 
independent of their arguments [65j . In this case, the Debye energy provides a natural cutoff for the integrals 
which enter the expression of the maximum Josephson current Jo (in zero voltage), namely, 

Jo = 2A L A R £ - F ]tk ( p ] \ > (143) 

kt ,kj?. 

and the integrals are evidently finite even in three dimensions when one approximates |tkj,kfil ~ \tk L k R \ ■ 
[Equivalently, within the above approximation one may transform the wave-vector integrals into energy 
integrals, whereby the density of state is taken at the Fermi level consistently with the BCS (weak-coupling) 
limit, in such a way that the ensuing integrals converge like in two dimensions.] 

It is, however, clear that this approximation does not apply when considering the BCS-BEC crossover, for 
which no upper momentum cutoff can be introduced in order to account for the formation of the composite 
bosons. In this case, the integral on the right-hand side of Eq. (|143p would be ultraviolet divergent, unless 
the factor |iki,k H | 2 therein vanishes rapidly enough for large |k/,| and |k^|. Making that integral to converge 
would thus require a resonable knowledge of the tunneling properties for fermions impinging on the barrier 
with large energies, knowledge which could only come from a separate study of the scattering on the barrier of 
superfluid fermions with all possible energies. This study is naturally provided by approaching the problem 
directly in terms of the BdG equations, as we have consistently done in this paper, thus making it useless 
any further recourse to the tunneling Hamiltonian approach. 
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